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POLYHEDRAL OMEGA: A NEW ALGORITHM FOR SOLVING LINEAR DIOPHANTINE SYSTEMS 


FELIX BREUER AND ZAEEIRAKIS ZAFEIRAKOPOULOS 


Abstract. Polyhedral Omega is a new algorithm for solving linear Diophantine systems (LDS), i.e., for 
computing a multivariate rational function representation of the set of all non-negative integer solutions to 
a system of linear equations and inequalities. Polyhedral Omega combines methods from partition anal¬ 
ysis with methods from polyhedral geometry. In particular, we combine MacMahon’s iterative approach 
based on the Omega operator and explicit formulas for its evaluation with geometric tools such as Brion 
decompositions and Barvinok’s short rational function representations. In this way, we connect two recent 
branches of research that have so far remained separate, unified by the concept of symbolic cones which 
we introduce. The resulting LDS solver Polyhedral Omega is significantly faster than previous solvers based 
on partition analysis and it is competitive with state-of-the-art LDS solvers based on geometric methods. 
Most importantly, this synthesis of ideas makes Polyhedral Omega the simplest algorithm for solving lin¬ 
ear Diophantine systems available to date. Moreover, we provide an illustrated geometric interpretation of 
partition analysis, with the aim of making ideas from both areas accessible to readers from a wide range of 
backgrounds. 
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1. Introduction 


In this article we present Polyhedral Omega, a new algorithm for computing a multivariate rational 
function expression for the set of all solutions to a linear Diophantine system. This algorithm connects 
two branches of research - partition analysis and polyhedral geometry - between which there has been 
little interaction in the past. To make this paper accessible to researchers from either field, as well as 
to readers with other backgrounds, we give an elementary presentation of the algorithm itself and we 
take care to motivate the key ideas behind it, in particular their geometric content. To begin, we use this 
introduction to define the problem of computing rational function solutions to linear Diophantine sys¬ 
tems and to give an overview of the algorithms developed in partition analysis and polyhedral geometry 
to solve it, before pointing out the benefits of Polyhedral Omega. 


Linear Diophantine Systems emd Rationed Functions. Let A e Z'”’'" be an integer matrix and let b e 
Z"* be an integer vector. In this article, we are interested in finding the set of non-negative integer 
vectors x e Z”q such that Ax S b. Since we are restricting our attention to non-negative solutions x e 
Z"q, it is equivalent to consider systems consisting of any combination of equations and inequalities. 
However, to streamline this article, we are going to focus on the case Ax S: b. We call such a system of 
constraints, given by A and b, a linear Diophantine system (LDS). 

Linear Diophantine systems are of great importance, both in practice and in theory. For example, 
the Integer Programming (IP) problem - which is about computing a solution to an LDS that maximizes 
a given linear functional - plays a pivotal role in operations research and combinatorial optimization 
[55]. In this article, we are not interested in finding one optimal solution, however. Instead, we wish to 
compute a rational function representation of the set of all solutions to an LDS. 

Of course, the set of solutions to a linear Diophantine system may be infinite. For example, the 
equation xi - X 2 = 0, which is equivalent to the system xi - X 2 ^ 0 and x\- X 2 ^ 0, has the solution 
set {(0,0), (1,1), (2,2),...}. One way to represent such infinite sets of solutions is via multivariate rational 
functions. Weidentify avector jc £ Z” with the monomial zp-z^^-.-.-z^" in n variables. Then, using 

the geometric series expansion formula, we can represent the above set of solutions by the rational 
function 
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1-2(1,1) I-Z 1 Z 2 

In general, we represent a set S c Z”q of non-negative integer vectors by the multivariate generating 
function 


(ps(z] = • 2 ''- 

xeS 


i.e., the coefficient z^ in (/>s is 1 if x e S and 0 if x ^ S. It is a well-known fact that when S is the set of 
solutions to a linear Diophantine system, then (ps is always a rational function ps- It is this rational 
function ps that we seek to compute when solving a given linear Diophantine system. We are not going 
to be interested in the normal form of ps, though, since the normal form may be unnecessarily large. 
Take for example the system Xi - 1 - X 2 = 100. Here, the set of solutions is 


S= {(100,0),(99,1),...,(0,100)} 


which can be represented by the rational function 

„(100,0) _(-l,101) 

f ^ ^ -I- -t 

l-2(-l,l) 

Note that the polynomial on the right-hand side is the normal form of this rational function, which 
has 101 terms. (In fact, the normal form of a rational function representing a finite set will always be 
a polynomial.) The rational function expression on the left-hand side is not in normal form since the 
denominator divides the numerator, however it is much shorter, having only 4 terms. We are therefore 
interested in computing multivariate rational function expressions, which are not uniquely determined, 
as opposed to the unique rational function in normal form. Given this terminology and notation, we can 
now state the computational problem of solving linear Diophantine systems that this article is about. 
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Problem 1.1. Rational Function Solution of Linear Diophantine System (rfsLDS) 

Input; A e Z"*’'", beZ’^ 

Output: An expression for the rational function p e Q (zi,..., z„) representing the set of all non-negative 
integer vectors x e Z"q such that Ax ^ b. 

Such rational function solutions to LDS are of great importance in many applications. For example, 
they can he used to prove theorems in number theory and combinatorics [4, 5, 50], compute volumes 
[16], count integer points in polyhedra ]17], to maximize non-linear functions over lattice points in 
polyhedral [32], to compute Pareto optima in multi-criteria optimization [30], to integrate and sum 
functions over polyhedra [12], to compute Grdbner bases of toric ideals [29], to perform various op¬ 
erations on rational functions and quasipolynomials [15], and to sample objects from polyhedra [52], 
from combinatorial families [35] and from statistical distributions [28]. We recommend the textbooks 
[13, 20, 31] for an introduction. 

Note that, while above rfsLDS is stated in terms of pure inequality systems, this restriction is not 
essential. The methods and algorithm we present in this article can be easily extended to handle the 
case of mixed systems directly, as does our implementation [23]. 

Partition Analysis. One of the major landmarks in the long history of Problem 1.1 is MacMahon’s semi¬ 
nal work Combinatory Analysis [50], published in 1915, where he presents partition analysis as a general 
framework for solving linear Diophantine systems in the above sense, particularly in the context of par¬ 
tition theory. The general approach of partition analysis is to employ a set of explicit analytic formulas 
for manipulating rational function expressions that can be iteratively applied to obtain a solution to a 
given linear Diophantine system. We will examine partition analysis in detail in Section 2; for now, we 
will confine ourselves to a quick preview to convey the general flavor. 

Consider the linear Diophantine inequality 

(1) 2xi-I-3x2 - 5^3 S’ 4 

To solve this system, MacMahon introduces the Omega operator defined on rational functions in 
IK(xi,X 2 ,X 3 )(A) in terms of their series expansion by 

+00 +00 

Os ^ 

5=-00 5=0 

where Cg e K(xi,X 2 ,X 3 ) for all s. In other words, acts by truncating all terms in the series expan¬ 
sion with negative A exponent and dropping the variable A. Using the Omega operator, the generating 
function fs of tho set of solutions S of (1) can be written as 

~ (1 - ZiA2)(l - Z2A3)(1 - Z3A-5) ■ 

The expression on the right-hand side of this equation is known as the crude generating function. 

Partition analysis is about developing calculi of explicit formulas for evaluating the Omega operator 
when applied to rational functions of a certain form. MacMahon’s 1st Rule provides a typical example: 

^11 
^^(l-Ax)(l-i) “ (l-x)(l-x*y)' 

Applied iteratively, such formulas allow us to symbolically transform the crude generating function into 
a rational function expression for (pgy desired. 

At the turn of the 20th century, Elliott was the first to develop such a calculus that yielded an algo¬ 
rithmic method for solving a linear Diophantine equation. At the same time, MacMahon was pursuing 
the ambitious and visionary project of applying such methods to partition theory. While algorithmic, 
Elliott’s method was not practical for MacMahon’s purposes, as it was too computationally intensive to 
carry out by hand. Therefore, MacMahon developed and continually extended his own calculus which 
enabled him to solve many partition theoretic problems, even though it did not constitute a general 
algorithm. Despite many successes, MacMahon did not achieve the goal he had originally set himself 
of using his partition analysis calculus to prove his famous formula for plane partitions [5]. It is perhaps 
for this reason that partition analysis fell out of interest for much of the 20th century, before it started 
to attract renewed attention in recent decades, beginning with Stanley [56]. 
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Figure 1. An overview of the field of partition analysis. 


It was Andrews who observed the computational potential of MacMahon’s partition analysis and 
waited for the right problem to apply it. In the late 1990s, the seminal paper of Bousquet-Melou and 
Eriksson [21] on lecture-hall partitions appeared. Andrews suggested to Paule to explore the capacity 
of partition analysis combined with symbolic computation. This collaboration gave rise to an ongo¬ 
ing series of over 10 papers and includes Omega [4] and Omega2 [5], two fully algorithmic versions of 
partition analysis powered by symbolic computation. This line of research has also attracted renewed 
interest in this field, including for example, a sequence of papers by Xin [64, 63, 65] who develops al¬ 
ternative partition analysis calculi and corresponding software packages. An illustration of the connec¬ 
tions between these different solutions to the problem is given in Figure 1. 

Polyhedral Geometry. Independently from the work on partition analysis, the field of polyhedral ge¬ 
ometry has made major contributions to the solution of linear Diophantine systems. Integer linear 
programming is a huge and very active field of research concerned with the computation not of a ra¬ 
tional function representation of the set of all non-negative integer vectors satisfying a given LDS, but 
instead with the computation of a single such vector which is optimal with respect to a linear functional. 
The integer programming problem is distinct from the one we are interested in in this article. Nonethe¬ 
less, two results on integer programming are important to mention here. On the one hand, deciding 
satisfiability of linear Diophantine systems is NP-complete, as many important and practical decision 
problems can be modeled easily using linear Diophantine systems [54, 55]. On the other hand, Lenstra 
has shown in 1983 [47] that if the number of variables is fixed, then the satisfiability of a linear Dio¬ 
phantine system can be decided and (if it is satisfiable) an optimal solution can be found in polynomial 
time. 

The problem of computing a rational function solution to a linear Diophantine system (rfsLDS), 
which we are concerned with in this article, is also at least as hard as any problem in NP, in the following 
sense. It is easy to see that NP-hard classes of the satisfiability problem for LDS can be reduced to the 
problem of computing the unique normal form rational function solution of the LDS. Interpolation ar¬ 
guments then show that moreover any rational function expression of polynomial size will readily yield 
the solution to the NP-hard problem in polynomial time. 

Thus, the question arises, whether rfsLDS also becomes polynomial time computable if the number 
of variables is fixed. Lenstra’s algorithm does not help here. The key obstacle is that a priori the encoding 
length of the output may be exponential in the encoding length of the input. As we have already seen in 
the example above, the normal form of the rational function solution of x -i- y = b has b+1 terms while 
the LDS itself has encoding length ^(logh). While in this particular case, we can find a rational function 
expression, —-, whose encoding length is in 0 (log b ), it is far from clear that such short rational 

function expressions exist in general. They do, however. 
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In 1993 Barvinok was able to show in a seminal paper [16] that for every LDS there does exist a ratio¬ 
nal function expression for its set of solutions whose encoding length is bounded polynomially in terms 
of the encoding length of the input - provided that the number of variables is fixed. Moreover, Barvinok 
gave an algorithm for computing such a short rational function representation. The key ingredient that 
makes these short decompositions possible is Barvinok’s method for computing a short signed decom¬ 
position of a simplicial cone into unimodular simplicial cones. These Barvinok decompositionswill play 
an important role in this article as well. 

Barvinok’s algorithm for solving rfsLDS combines Barvinok decompositions with a number of other 
sophisticated tools from polyhedral geometry, ranging from the computation of vertices and edge di¬ 
rections of polyhedra [40, 41] to the triangulation of polyhedral cones [34, 53]. This combination makes 
Barvinok’s algorithm very complex and much more involved than the relatively straightforward collec¬ 
tion of rules based on explicit formulas that MacMahon had envisioned in the partition analysis frame¬ 
work. A testament of this difficulty is the fact that even though Barvinok’s algorithm quickly grabbed the 
attention of scientific community, it took 10 years and a sequence of additional research papers, before 
Barvinok’s algorithm was first implemented in 2004 by De Loera et al. [33]. 

The importance of Barvinok’s short rational function representations in this area cannot be over¬ 
stated, as they have given rise to a whole family of algorithms and theoretical complexity results for a 
large range of problems, many of which we have mentioned above. Of particular note in our context is 
that Barvinok and Woods proved a theorem in [15] which implies in particular that MacMahon’s Omega 
operator can be evaluated in polynomial time in fixed dimension - not only on crude generating func¬ 
tions as they arise in partition analysis but on a wider class of rational functions. 

To conclude our overview of prior literature, we briefly mention a couple of further related publica¬ 
tions [25, 43, 44], a detailed discussion of which is out of scope of this article. [43, 44] give two different 
algorithms for solving rfsLDS via an explicit formula, phrased in terms of a summation over all sub¬ 
matrices of the original system. This can be viewed as a brute force iteration over all simplicial cones 
formed by the hyperplane arrangement given by the system, together with a recursive algorithm for 
computing the numerators of the corresponding rational functions. In contrast to an explicit formula 
given in [25], the formulas from [43, 44] are algorithmically tractable, though the running time can be 
exponential, even if the dimension is fixed. Since the iterative application of is not the focus of 
[43, 44], we do not view these algorithms as part of the partition analysis family for the purposes of this 
article. To our knowledge these algorithms have not been implemented. 


Our Contribution. In this article, we bring partition analysis and polyhedral geometry together. Some¬ 
what surprisingly, these two branches of research have so far had relatively little interaction. To bridge 
this gap, we provide a detailed and illustrated study of these two fields from a unified point of view. 
In particular, we interpret the partition analysis calculi provided by Elliott, MacMahon and Andrews- 
Paule-Riese from a geometric perspective. For example, MacMahon’s ansatz of setting up the crude 
generating function corresponds in the language of geometry to a clever way of writing an arbitrary 
polyhedron as the intersection of a simplicial cone with a collection of coordinate half-spaces - a con¬ 
struction which we call the MacMahon lifting. While the connection between rational functions and 
polyhedral cones is well-known and a discussion of the Elliott-MacMahon partition analysis algorithm 
from this point of view can he found in [56], some of this material is new, including our geometric dis¬ 
cussion of Andrews-Paule-Riese and Xin. 

The main contribution of this article is Polyhedral Omega, a new algorithm for solving linear Dio- 
phantine systems which is a synthesis of key ideas from both partition analysis and polyhedral geom¬ 
etry. Just like partition analysis algorithms. Polyhedral Omega is an iterative algorithm given in terms 
of a few simple explicit formulas. The key difference is that the formulas are not given in terms of ra¬ 
tional functions, but instead in terms of symbolic cones. These are symbolic expressions for simplicial 
polyhedral cones, given in terms of their generators. Applying ideas like Brion’s theorem from poly¬ 
hedral geometry gives rise to an explicit calculus of rules for intersecting polyhedra with coordinate 
half-spaces, i.e., for applying the Omega operator on the level of symbolic cones. 

After the Omega operator has been applied, the final symbolic cone expression can then be con¬ 
verted to a rational function expression using either an explicit formula based on the Smith Normal 
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Form of a matrix, or using Barvinok decompositions. The Smith Normal Form approach has the advan¬ 
tage of being simpler and easier to implement, while the Barvinok decomposition approach is much 
faster on many classes of LDSs and produces shorter rational function expressions. However, for many 
applications it can be advantageous to not convert the symbolic cone expression to a rational function 
expression at all, since this conversion can increase the size of the output dramatically. Especially when 
the output is intended for inspection by humans, we strongly recommend working with the symbolic 
cone expression as it is typically much more intelligible. 

Polyhedral Omega has several advantages. In comparison with partition analysis methods. Polyhe¬ 
dral Omega is much faster. In fact, on certain classes of examples it offers exponential speedups over 
previous partition analysis algorithms - even if Barvinok decompositions are not used! - by virtue of 
working with symbolic cones instead of rational functions and through an improved choice of transfor¬ 
mation rules. In comparison with polyhedral geometry methods. Polyhedral Omega is much simpler. It 
does not require sophisticated algorithms for explicitly computing all vertices and edge-directions of a 
polyhedron or for triangulating non-simplicial cones. These computations happen implicitly through 
the straightforward application of a few simple rules given by explicit formulas. At the same time, if 
Barvinok decompositions are used for the conversion of symbolic cones into rational functions, then 
Polyhedral Omega runs in polynomial time in fixed dimension, i.e., it lies in the same complexity class 
as the fastest-known algorithms for the rfsLDS problem. Polyhedral Omega is the first algorithm from 
the partition analysis family that achieves polynomial run time in fixed dimension. Moreover, Polyhe¬ 
dral Omega is the first partition analysis algorithm for which a detailed complexity analysis is available. 
This analysis is made possible by our geometric point of view. 

On the whole. Polyhedral Omega is the overall simplest algorithm for solving linear Diophantine 
systems, while at the same time its performance is competitive with the current state-of-the-art. 

In this article, we give a theoretical description and analysis of Polyhedral Omega, together with a 
detailed motivation and illustration of the central ideas that makes the exposition accessible to readers 
without prior experience in either partition analysis or polyhedral geometry. The definition of the al¬ 
gorithm is explicit enough to make it possible for interested readers to implement the algorithm even 
without being experts in the field. Moreover, we have implemented Polyhedral Omega on top of the 
Sage system [57] and the code is publicly available [23]. 

Organization of this Article. In Section 2, we give a geometric interpretation of previous partition anal¬ 
ysis calculi by Elliott, MacMahon, Andrews-Paule-Riese and Xin. In Section 3, we present the Polyhedral 
Omega algorithm and prove its correctness. We take particular care to motivate and illustrate the key 
ideas in an accessible manner. In Section 4, we give a detailed complexity analysis of our algorithm. 
In Section 5, we then compare Polyhedral Omega to other algorithms from both partition analysis and 
polyhedral geometry. Finally, we conclude in Section 6, where we also point out directions for future 
research. 


2. Partition Analysis and its Geometric Interpretation 

In this section we present some of the major landmarks of partition analysis and interpret the results, 
which are typically phrased in the language of analysis, from the perspective of polyhedral geometry. 

In the following, we will introduce the required concepts and terminology from polyhedral geometry 
as we go along. However, a comprehensive introduction to polyhedral geometry is out of scope of this 
article. We therefore refer the reader to the excellent textbooks [20, 31, 54, 66] for further details. 

2.1. MacMahon. 

The Qj operator and the crude generating function. MacMahon presented his investigations concern¬ 
ing integer partition theory in a series of (seven) memoirs, published between 1895 and 1916. In the 
second memoir [49], MacMahon observes that the theory of partitions of numbers develops in parallel 
to that of linear Diophantine equations. In particular, he defines the operator in Article 66: 

“Suppose we have a function P{x, a) which can be expanded in ascending powers of x. Such ex¬ 
pansion being either finite or infinite, the coefficients of the various powers of x are functions of a 
which in general involve both positive and negative powers of a. We may reject all terms contain¬ 
ing negative powers of a and subsequently put a equal to unity. We thus arrive at a function of x 
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only, which may be represented after Cayley (modified by the association with the symbol by 
F{x, a) the symbol S denoting that the terms retained are those in which the power of a is S 0.” 
A modern wording of the definition is given in [4]. The operator is defined on functions with abso¬ 
lutely convergent multisum expansions 

OO OO OO 

Y. Y Y '"^r 

Si=—ooS2=—oo Sr=—oo 

in an open neighborhood of the complex circles \Xi\ - 1. The action ofD.^ is given by 

OO OO OO OO OO OO 

Y Y ■■■ Y Y Y 

Si=-00S2=-00 Sr=-00 Si=0S2=0 5^=0 

MacMahon proceeds with defining what he calls the crude generating function. This is a (multivariate) 
generalization of an idea appearing in Cayley as well as in Elliott’s work (see next section). Given a linear 
system of Diophantine inequalities, MacMahon introduces an extra variable for each inequality. These 
extra variables are denoted by A. Given A e and b e Z™ we consider the generating function 

m 

<,,(z,A)= x 


and based on the geometric series expansion formula (1 - z) ^ we can transform the series 

into a rational function. The rational form of <1)^ j^fz) is denoted by ^(z), i.e.. 


p^_^(z,A) = A *n 


/=! fl-ZiA(^^)^' ' 
. 


We call p^^(z, A) the X-generating function of A and b. The crude generating function is the syntactic 
expression 




A-^n 


i=ifl-ZiA(^^)i' 

. 


obtained by prepending to the A-generating function. We present an example that we also use for 
the geometric interpretation. Let A = [ 2 3 ] and b - 5 and consider the Diophantine system Ax ^ b. 

Then 

(b\.{zi,Z 2 ,X)- y and p^ . (zi,Z 2 , A) =- 

' (1-ZiA2)(1-Z2A3) 


Thus, the crude generating function for this system is 

A-s 

(1-ZiA2)(1-Z2A3) 

After defining the crude generating function, MacMahon proceeds to evaluate this expression by the 
use of formal rules, such as 

1 _ 1 

a-Xx] {!-{,) ~ (l-M(l-x^y)' 

In what follows we will explain what this means geometrically. 


Cones and fundamental parallelepipeds. The connection between partition analysis and geometry hinges 
on the fact that certain rational function expressions correspond directly to generating functions of lat¬ 
tice points in polyhedral cones. A lattice point is any integer linear combination of a fixed set of basis 
vectors. Without loss of generality, in this article, we will always work with the the standard unit vectors 
as our basis, so that “lattice point” is synonymous with “integer point”. 

We already used O to denote the formal power series generating function for linear Diophantine sys¬ 
tems. More generally, for a set S c of lattice points we will use <l>s to denote the generating function 
of S, i.e., <l>s = Uses ■2*- For convenience, if S c is a set of real vectors, we will also write Cis to denote 
the generating function of all lattice points in S, namely <l>s = Hs^SnZ'^ Of course, if and where Os 
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(A) ((1,0), (1,3)). 
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(B) n ((1,0), (1,3)) decom¬ 
posed into three translates of 



(C) The corresponding tiling of 
((1,0), (1,3)) hy translates of its 
fundamental parallelepiped. 


Figure 2. Generating functions of cones. 


converges depends on the set S, but these questions will be of no particular concern to us, as we will see 
below. We are primarily interested in generating functions <I>c where C is a simplicial polyhedral cone. 
We define the cone C generated by vi,V 2 ,...,v^ as 

id 1 






0 ^ A; e IR L 


11=1 J 

The Vi are the generators of C. If the are linearly independent, the cone C is simplicial. All cones in 
this article are rational polyhedra, which means that they can also be defined as the set of real vectors 
satisfying a finite system of linear inequalities with integer coefficients. A cone C is pointed or line-free 
if it does not contain a line. All simplicial cones are pointed. For pointed cones C the power series Tic 
always has a non-trivial radius of convergence, see also Section 2.4. If a given cone C is pointed and 
the chosen set of generators Vi is inclusion-minimal, then the sets fi are called the extreme rays of 
the cone. The set of extreme rays i?; of a pointed cone is uniquely determined, but within each ray the 
choice of generators 0^ vte Ri is arbitrary. A vector v e Z” is primitive, if the greatest common divisor 
of its entries is 1. For every vector Of v e <Q" there exists a unique primitive vector primfy) that is a 
positive multiple of v. Equivalently, prim(i;) is the unique shortest non-zero integer vector on the ray 
V- With these normalizations we find that each pointed rational cone has a unique minimal set of 
primitive generators. 

Later, we are also going to need the notion of affine cones, which just means cones like those de¬ 
fined above translated by a rational vector q. For translates of pointed cones, the translation vector q is 
uniquely determined, in which case q is called the apex of the affine cone. We wUl frequently drop the 
adjective “affine” and just refer to cones in general. It wUl be clear from context whether a given cone 
has its apex at the origin or not. For a given q e we adopt the notation 


'ig^{vi,...,vd;q) ^q-i-'^^{vi,...,Vd). 


Flere, -i- denotes the Minkowski sum, i.e., for two sets A and B we have A-t B - {a-t b \ ae A,beB]. 

If instead of considering all conic combinations of the generators with real coefficients, we consider 
only non-negative integer combinations, then we obtain the semigroup 

1 d 1 


‘^^{Vi,...,Vd) 




0 ^ A; E Z L 


! = 1 


which we also call the discrete cone generated by the vi. As above we write ■- + 

^^(Vi,...,Vd). 

As an example, consider the semigroup generated by (1,0) and (1,3) as shown in Figure 2a. If we take 
the geometric series and expand it, i.e., 1 -i- z, -i- Zj -i-..., we observe that this series is the generating 
function of the lattice points on the non-negative x-axis. In order to take the generating function for 
the lattice points at height 3, we need to multiply the series by (ziz^). Accordingly, to take the lattice 


















POLYHEDRAL OMEGA: SOLVING LINEAR DIOPHANTINE SYSTEMS 


9 


points at height 6 we multiply hy and so on. It is easy to see, that the generating function of the 

semigroup is the product of two geometric series, namely and < so that we obtain 

1 

Following the same idea, we can translate the whole semigroup by multiplying by a monomial. If we 
multiply --jT by Z 1 Z 2 we obtain the generating function of the red points in Figure 2b. Similarly, 

we multiply by z\z\ to obtain the generating function of the green points. Now, since there are no 
overlaps, we can sum the generating functions in order to obtain the generating function for all lattice 
points in the cone generated by (1,0) and (1,3) with real coefficients, 

2'! 1 + Z 1 Z 2 +-Zl^l 

The numerator of this rational function encodes the lattice points (0,0), (1,1) and (1,2) which are pre¬ 
cisely the lattice points in the parallelogram 11 = {Ai(l,0) - 1 - A(l,3) | 0 ^ A; < 1} as shown in blue in Fig¬ 
ure 2c. Moreover, the cone ((1,0), (1,3)) is tiled by translates of 11. This observation generalizes. 
Given a simplicial cone C = {vi,V 2 ,..., Vd) e IPS'*, we define its fundamental parallelepiped as 

[ d 

n„(C) = n„(yi,...,yrf) = 0^ A/< 1,A, eR 

and the set oflattice points in the fundamental parallelepiped as n^^(C) = YL^{C)nZ‘^. Given this defini¬ 
tion, the general form of the generating function of a simplicial cone C . .,Vm) generated 

by fi,.... I'm c is 


( 2 ) 


<I'f.{Zi,...,Zd) 


Z“ 


where ^^^[zi,...,Zd) is the generating function for the lattice points in the fundamental paral¬ 
lelepiped of C. This fundamental observation goes back to Ehrhart [36, 37, 38]. From this identity it 
is immediately obvious that cones C with few lattice points in their fundamental parallelepiped 11^^ (C) 
have particularly short representations as rational functions. Cones C where 11^^ (C) contains just a 
single lattice point are called unimodular. 

An important property of the fundamental parallelepiped is that its lattice points are coset represen¬ 
tatives for the lattice points in the cone modulo the semigroup of the generators (as exemplified by the 
construction of the generating function). To put it in another way: the set of lattice points in the cone 
is the set of lattice points in the fundamental parallelepiped, translated by all non-negative integral 
combinations of the generators of the cone, i.e.. 


(3) GHZ'* = . . U L*>; + n2.(C) = '^2(yi,...,y„) + n^,(C) 

Here U denotes an ordinary union together with the assertion that the operands are disjoint. 


Geometry of the operator and theMacMahon lifting. Now we can see the connection of MacMahon’s 

construction to polyhedral geometry. In Figure 3, we see the geometric steps for creating the A-cone 
(equivalent of the A-generating function) by means of the MacMahon lifting and then applying the 
operator. 

Given an LDS Ax ^ b with m inequalities and d variables, we first lift the standard generators of 

the positive orthant of by appending the exponents of the A variables (thus lifting the generators in 

This lifting idea first appears in an analytic context in Cayley’s work [27]. Consider the cone 

generated by the MacMahon lifting, i.e., the cone generated by the columns of matrix V - ('jf j- The 

generating function of this cone, due to (2), is exactly - ttitt- This is true because the matrix V is 

jl-Zil'-'’' h'J 

unimodular by construction, whence the fundamental parallelepiped of the cone contains exacdy one 
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(A) Consider the first 
quadrant and lift the 
standard generators 
according to the input. 



(b) Shift according to 
the inhomogeneous 
part in he negative 
A direction the cone 
generated by the lifted 
generators. Intersect 
with the xy-plane. 



(C) Take the part of 
the cone (a polyhe¬ 
dron) living above the 
xy-plane and project 
its lattice points. 


A 



(d) The projected 
polyhedron contains 
the solutions to the 
inequality. 


Figure 3. as lifting and projection. 


lattice point, namely the origin. Next, we translate the cone by q - (_°^) and denote the new cone by 
C = ^(V) -t- Then the generating function of cone C is exacdy the A-generating function, i.e., 

m 2 

Applying the Omega operator to <I>c now corresponds geometrically to intersecting C with all half¬ 
spaces where the A-variables are non-negative and then applying a projection that forgets the A-variables. 
To make this precise, let us denote by Hi the intersection of coordinate halfspaces where all A variables 
are non-negative, i.e., 

Ha = |(xi, X2,...,Xd, xaj , xaj, ..., Xij e | x^. S 0 for 1 ^ ^ mj-. 

Next, let 71 denote the projection from to that forgets the last m coordinates which correspond 
to A-variables. Given this notation, the action of on the A-generating function corresponds to first 
intersecting C with the Hi and then projecting the resulting polyhedron using 7i, i.e., 

Continuing with the example we had chosen previously, assume A = [2 3] and h = 5. In Figure 3a, we 
construct the vectors (1,0,2) and (0,1,3), lifting the two standard basis vectors (1,0) and (0,1) according 
to the entries of matrix A. In Figure 3b, we shift the cone generated by (1,0,2) and (0,1,3) by -b = -5, 
and consider the halfspace defined by the xy-plane, where A is non-negative. The two intersection 
points of the cone with the hyperplane are marked. Now, in Figure 3c, we consider the intersection of 
the halfspace defined above with the shifted cone. The lattice points in the intersection are shown in 
blue marked, and their projections onto the A = 0 plane are shown in red. These projections are the 
non-negative solutions to the original inequality, as shown in Figure 3d. 

MacMahon’s rules. MacMahon’s method was general and in principle algorithmic, based on Elliott’s 
work (see next section). As we shall see, Elliott’s method is simple but computationally very ineffi¬ 
cient. MacMahon introduced a set of rules in order to deal with particular combinatorial problems 
computationally - in particular, as he had to carry out his computations by hand. We present some of 
MacMahon’s rules from a geometric point of view. 
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Figure 4. The geometry of the first MacMahon rule. 



Figure 5. The geometry of the second MacMahon rule. 


The first two rules MacMahon gave are 

1 


and 




n> 


(1-Ax)(l-^) (l-x)(l-x*y) 


1 , ^ 


for s G ] 


for s G; 


->o 


-?o* 


(1-A^x)(l-|) (l-x)(l-xy*) 

There is an apparent symmetry in the input, hut elimination results in structurally different numerators, 
1 — 

i.e., 1 versus 1 + xy . This can he easily explained hy recalling that the numerator of the generat¬ 
ing function enumerates the lattice points in the fundamental parallelepiped. In Figures 4 and 5 the 
respective fundamental parallelepipeds are shown. 

The fourth rule in MacMahon’s list is 

1 


o> 


1-xyz 


(1 - Ax) (1 - Ay) (1 - I) (1 - X) (1 - y) (1 - xz) (l - yz)' 


Although the rational function on which acts has three factors in the denominator, the resulting 
rational generating function has four factors in the denominator and the numerator has terms with 
hoth positive and negative signs. This is because the initial cone defined by the inequality x-i- y - z 0 is 
not simplicial. The cone admits a signed decomposition into simplicial cones. When summing up the 
generating functions of the cones we obtain the right hand side of the rule, with a negative term due to 
the signed decomposition. The decomposition is shown in Figure 6. 
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Figure 6. The fourth MacMahon rule. We start with the non-simplicial 
cone ((1,0,0), (0,1,0), (1,0,1), (0,1,1)). Then we consider the simplicial cones 
((1,0,0), (0,1,0), (1,0,1)) andS^jj ((0,1,0), (1,0,1), (0,1,1)) as well as their intersection 
SfjjdO, 1,0),(1,0,1)). By inclusion-exclusion we have a decomposition of the initial 
cone. 


2.2. Elliott. One of the first references relevant to partition analysis is Elliott’s article “On linear homo¬ 
geneous Diophantine equations” [39]. The work of Elliott is exciting, if not for anything else, because 
it addresses mathematicians that lived a century apart. It was of interest to MacMahon, who based his 
method on Elliott’s decomposition, but also to 21st century mathematicians for explicitly giving an im¬ 
portant algorithm. Although the algorithm has very bad computational complexity, it can be considered 
as an early algorithm for lattice point enumeration (among other things). 

The problem Elliott considers is to find all non-negative integer solutions of one homogeneous linear 
Diophantine equation 

m m+n 

(4) ^ UiXi - ^ biXi = 0 for a,-, h,- e Z^o- 

i=l i=m+l 

The starting point for Elliott’s work is the fact that even if one computes the set of “simple solutions”, 
i.e., solutions that are not combinations of others, there are syzygies preventing us from writing down 
formulas giving each and every solution to the equation exactly once. He proceeds by explaining that 
his method computes a generating function whose terms are in one-to-one correspondence with the 
solutions of the equation. 

The basic idea employed by EUiott (and the basis of partition analysis) is the introduction of an extra 
variable denoted by A, encoding the equation. This idea, as we saw, was employed by MacMahon and 
can be traced back to Cayley. Elliott then observes that the terms in the series expansion of 

1 

(1 - XlA“l) • (1 - (l - Xm+lA-^fn+l) ■ (l - XlA-^'n+n) 

correspond to solutions of (4). The question then is how to decompose this rational function into a 
sum of rational functions, such that the expansion of each of them contains either only terms with 
non-negative A exponents or only terms with non-positive A exponents. Given such a decomposition 
it is safe to drop the components that give terms containing A. Elliott’s algorithm computes such a 
partial fraction decomposition. Given a rational function expression Of where m; is a monomial 
Zi, Z 2 ,..., Zfc, A, A“^ and k e Z^o, Elliott computes a sum of the form X + Z with pi being 

a monomial in zi,Z 2 ,.. .,Zfc, A and qj a monomial in zi,Z 2 ,.. .,Zfc, A”'^, based on the identity 

1 1 11 

(5) - =-1- 

(l-ziA“)(l-||) (l-ziZ 2 A“-/*)(l-ziA“) (l-ziZ 2 A“-/*)(l-||) [I- 

To achieve the goal of only having summands where the exponents of A in the denominator all have 
the same sign, the above identity has to be applied iteratively. The left-hand side of (5), with exponents 
[a,-p), is reduced to (a - p, a) or (a - p,-p), depending on whether p. This leads to a recurrence 
corresponding exactly to the Euclidean algorithm, see also [22]. 

Elliott’s identity has a very nice geometric interpretation, which is illustrated in Figure 7. We observe 

that-^ is the generating function of the 2-dimensional cone C = f(l,0, a), (1,0,-d)1, 

(l-xA“)(l-^) f RV 

i.e., the fundamental parallelepiped of C contains only the origin. Since the point v = (1,1, o: - /3) is 
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(C) Drop green cone 



(D) ((1,0,1), (1,1,-2)), 

( 2 , 1 ,- 1 ) 



Figure 7. Elliott’s decomposition method. 


in the interior of the cone C, the cones A - ((1,0, a], (1,1, a - /!)) and B - ((1,0,-/I), (1, 1, a - /!)) 

subdivide cone C. Their intersection is exactly the ray starting from the origin and passing through 
V - {1,1, a - P). By a simple inclusion-exclusion argument we have the signed decomposition C = A + 
B - {An B). This decomposition translated hack to the generating functions level is exactly the partial 
fraction decomposition employed hy Elliott. One should mention here that Elliott insists on the fact 
that the numerators are +1, which geometrically means that he deals with unimodular cones, which 
also play a key role in Barvinok’s algorithm. 

2.3. Andrews-Paule-Riese. In 1997, Bousquet-Melou and Eriksson presented a theorem on lecture hall 
partitions in [21]. This theorem gathered a lot of attention from the community (and it still does, with 
many lecture hall type theorems appearing still today [18]). Andrews, who had already studied MacMa- 
hon’s method and was aware of its computational potential, figured that lecture hall partitions offered 
the right problems to attack algorithmically via partition analysis. 

At the same time it was planned to spend a semester during his sabbatical at RISC in Austria to work 
with Paule. It is only natural that the result was a fully algorithmic version of MacMahon’s method 
powered by symbolic computation. This collaboration gave a series of 10 papers [4, 1, 10, 5, 6, 7, 3, 8, 9, 
2]. Many interesting theorems and different kinds of partitions are defined and explored in this series of 
papers, but for us the most important two are [4] and [5], which contain the algorithmic improvements 
on partition analysis. Namely, in [4] the authors introduce Omega, a Mathematica package based 
on a fully algorithmic version of partition analysis, while in [5] they present a more advanced partial 
fraction decomposition (and the related Mathematica package Omega2) solving some of the problems 
appearing in Omega. While presenting the methods, we will see some of their geometric aspects. 

The main tool in Omega is the Fundamental Recurrence (Theorem 2.1 in [4]) for the operator. 
Following MacMahon, or Elliott for that matter, iterative application of this recurrence is enough for 
computing the action of O^. The interested reader can find the definition of the recurrence in [4]. The 
fundamental recurrence assumes that the exponents of A in the denominator are +1. This is not a strong 
assumption, as noted in [4], since we can always employ a decomposition to linear factors. The obvious 
drawback of this approach is that we introduce complex coefficients instead of +1. This motivates the 
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(A) Assignment of index subsets Ip to faces F. The 
ray generated by is indexed by / = {1,2}. The face 
generated by fli, 02 is indexed by / = {2). 


y 



(b) The simplicial cone {ai,a 2 ,a 3 ) has one 

open facet, namely the red facet on which the coef¬ 
ficient of 02 is zero. The vector (0,1,0) denotes that 
the facet with index set I = {2} is open. Similarly, 
the cone [ 01 , 02 , 03 ] has both the hlue and 

red facets open, which correspond to the index sets 
{1} and {2}, respectively. 


Figure 8. Half-open cones. 


desire for a better recurrence. In [5] the same authors introduce an improved partial fraction decompo¬ 
sition method (the generalized partial fraction decomposition), given by the following recurrence.'^ Let 
gcd(a, ;S) = 1 and inv^(a) denote the multiplicative inverse of a modulo b. Then 


( 6 ) 


where P, 


1 


a.p 


Qa,j6 


a ,/6 • 


(l-ZiZ“)(l-Z2-Z3) (1 

a-1 )S-1 

^ a,Z 3 and ;= ^ hiZg for 


i=o 


!=0 


■Zj^z“)(l- 


Z 2 Z 3 ) 


(7) 


(inv^(a)i mod (inv„(/i)i mod a) 

Z, Zr, 


( 8 ) 


bj - 


if p\i or i - 0, 
otherwise, 
if i = 0, 
otherwise. 


1 "2 

„a 

1 ^2 

(inv^(a)i mod (inv^yl)/ mod a) 

1^1 ^2 

Before interpreting the rational functions involved in the Fundamental Recurrence of [5] as cones in 
Figure 9, we introduce some notation and the concept of a half-open cone, which is useful when dealing 
with cone decompositions. The dimension of a polyhedron is the dimension of its affine hull. A face of 
a polyhedron P is a subset of P on which some linear functional is maximized. Faces of polyhedra are 
polyhedra themselves. Faces capture the intuitive notion of a “side” of a polyhedron. 0-dimensional 
faces are called vertices, 1-dimensional faces are called edges and the [d- l)-dimensional faces of a 
polyhedron of dimension d are called facets. For example, a 3-dimensional cube has 8 vertices, 12 
edges and 6 facets. Recall that 


C^^^{vi,V2,...,Vk;q)^ q+\ J^i\-iVi 

i=l 


0 ^ A E [ 


Assume C is simplicial and let [A:] = 7 c Z^o be the index set for its generators. Then each face F of C 
can be identified by a set Ip c /, since a face has the form 


F - q + 


k 

E 

1=1 


0 ^ A; e IR, Xj - 0 for j e Ip 


Essentially, this means that the face is generated by a subset of the generators. Note that Ip contains the 
indices of the generators of the simplicial cone C that are not used to generate F as in Figure 8a. A facet 
is identified with the index of the single generator not used to generate it. 


'Here, we have slightly rewritten the recurrence from [5| to facilitate the geometric interpretation. 
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(C) Cone A 


(d) Cone B 


(e) The decomposition 


Figure 9. Omega 2 decomposition. 


A half-open cone C is a cone from which some facets have heen removed. From the previous discus¬ 
sion, for each facet we need only to mention the generator corresponding to that facet. We will use a 
0-1 vector o to express the openness of a cone, i.e., the openness vector has an entry of 1 in the position 
k if the facet corresponding to the fc-th generator is open. To this end, we define 




d 

i=l 


0 ^ Aj e IR if Oi-Q and 0 < A,- e IR if ot-l 


The above notation will he used throughout the paper. To improve readability, we will drop the o from 
the notation when o = 0, i.e., when all facets are closed. Similarly, we drop q when q is the zero vector. 
Moreover, if V is a matrix with vectors ,..., as columns we may write V instead of listing the vectors 
explicitly, such that, for example, {V ),..., 0). 

Note that both (2) and (3) extend naturally to the case of affine half-open cones. Let V be the matrix 
with fi,..., yrf E Z" as columns. If we generalize the notion of a fundamental parallelepiped to 




0 ^ A/ < 1 if Oi-0 and 0 < A,- ^ 1 if o,=l 


and let n°„{V) := n“(l^) n Z”, then we have both 


(9) 






ntid- 




and <i^°{V;q)nZ"^q^^{V;q) + nl„{V). 


Returning to the cone interpretation of Omega2 (Figure 9), we observe that the 2-dimensional cone 
C = Sg"jj((l,0,a),(0,1,/i)) is the left hand side of the generalized partial fraction decomposition (6). 
Then we intersect the xy-plane with the plane the cone lives in. Choose a ray on this intersection 
(conveniently, one such ray is {-f, a,0), denoted by red in the figure). Then consider two cones A - 
a,0),(l,0,a)) and B - ((-j0, a,0), (0,1,/i)). We observe that the cone C is decomposed 

as the cone A minus the half-open cone B, as can be seen in Figure 9. In particular, the generating 
function for the fundamental parallelepiped of cone A is exactly the numerator and similarly the 
fundamental parallelepiped of B is precisely the numerator Qa,p- See Figure 10 for an example. This 
connection between polynomials of this form and fundamental parallelepipeds is weU-known in the 
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Figure 10. Rule (6), for a = 13 and p - 5, decomposes C = ‘i^(gi,g 2 ) into the differ¬ 
ence of = '^^(g,gi) andS = i-e., [C] = [Al] - [R], where gi = (1,0,13),g2 = 

(0, l,5),g = (-5,13,0). While the fundamental parallelepiped of C contains just a sin¬ 
gle lattice point, the lattice points in the fundamental parallelepipeds of A, B corre¬ 
spond to the numerators Pa,p and Qa,p respectively. The number of lattice points is 
linear in a, p, i.e., it is exponential in the encoding length of a, p. 


context of Dedekind-Carlitz polynomials [19]. In particular, these polynomials have a short recursive 
description [22], which can be exponentially more compact than (7) and (8). 

This shows that the Andrews-Paule-Riese decomposition adds a layer of complexity to the geometry 
of partition analysis in comparison with previous methods. Elliott is dealing only with unimodular 
cones, but has to resort to a recursive procedure for the elimination of a A variable. In contrast, the 
generalized partial fraction decomposition of [5], eliminates a A variable in one step at the expense 
of having cones with fundamental parallelepipeds containing potentially exponentially many lattice 
points. Moreover, the Andrews-Paule-Riese decomposition is more complex than previous rules in that 
it works with half-open cones. 


2.4. Questions of Convergence. In the partition analysis series of papers by Andrews-Paule-Riese and 
in particular in [4], the authors stress the importance of working with absolutely convergent series. This 
is because the Omega operator is not well defined if rational function expansions do not have a common 
region of convergence. Quoting from [4]: 

“While MacMahon did not carefully distinguish whether his Laurent series were ana¬ 
lytic or merely formal, we emphasize that it is essential to treat everything analytically 
rather than formally because the method relies on unique Laurent series representa¬ 
tions of rational functions. For instance, if we were to proceed formally, then 


n> 


n=C\ n=0 ^ 


while - Q; 


n=0 


■■ E 

n=l 


’A-" = 0. 


But if we allowed a purely formal application of the geometric series, then both initial 
expressions are .” 

While this observation about the Omega operator is indeed crucial, the conclusion that it is therefore 
“essential to treat everything analytically rather than formally” is not valid, as there are several methods 
to treat these issues from a formal point of view. 

Xin, who also gives an improved set of rules for computing the Omega operators in [64], notes that 
by working in the field of iterated Laurent series and by normalizing the rational function expression 
accordingly, it is possible to disregard questions of convergence. An in-depth study of different ways to 
define multivariate Laurent series appears in [11]. There, Aparicio Monforte and Kauers introduce mul¬ 
tivariate Laurent series as linear combinations of series that are supported in possibly shifted pointed 
cones compatible with a given ordering. They prove that this construction ensures that the support is 
well ordered and embed their theory in the Malcev-Neumann series construction presented by Xin in 
[64]. Their construction is a natural and “large” field of multivariate Laurent series. 
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Independently from these developments, researchers in the polyhedral geometry community have 
heen dealing with these analytical phenomena from a geometric perspective for a long time. Recall that, 
intuitively, we think of as corresponding to the I-dimensional cone [0,+oo) and of as corre¬ 
sponding to the half-open i-dimensional cone (- 00 , 0). The identity on the level of rational 

functions then leads us to considering equivalence classes of cones “modulo lines”, see Section 3. This 
is a routine practice in the polyhedral geometry literature [13] which is used for great effect (for example 
in Brion’s famous theorem, see Theorem 3.1 below) and should therefore be regarded as “a feature not 
a bug”. 

In this article, we have to be careful, however, since, as the above example shows, the Omega operator 
is not well-defined on equivalence classes modulo lines. Nonetheless, as we will see in the next section, 
polyhedral methods can still be applied to evaluate the Omega operator in an entirely symbolic manner, 
as long as we work with the right representatives of these equivalence classes modulo lines. To choose 
representatives consistently, we define any cone C as forward if all generators of C are lexicographically 
larger than the origin, i.e., if their first non-zero entry is positive. The notion of forward cones appears in 
[45] and plays a prominent role in the Lawrence-Varchenko decomposition (Theorem 3.2) that we use 
below. For the Omega operator to be well defined, we require that all series appearing in the expression 
on which Omega acts are supported in the half-open half-space of vectors that are lexicographically 
larger or equal than the origin. In general, one can define forward with respect to an arbitrary half-open 
half-space, but the above definition will suffice for our purposes. By working with expansions of rational 
functions that are forward in this geometric sense, we can apply the Omega operator without recourse 
to a full analytic treatment. In the example given above, the series is not a forward ex¬ 

pansion of Note that Xin’s iterated Laurent series field construction corresponds precisely to this 
geometric idea of working with linear combinations of cones that are all forward in the lexicographic 
sense, again see [11]. In Section 3, we will discuss both the method of working with cones modulo lines 
as well as the notion of forward cones in detail. 

2.5. Applying Omega vs. Solving LDS. MacMahon invented the operator for solving linear Dio- 
phantine systems. In this framework, is applied to rational functions that arise from the MacMahon 
lifting. However, once is defined, the question arises how to evaluate when applied to a general 
rational function r whose denominator factors into binomials. It turns out that this seemingly more 
general problem can be reduced to solving linear Diophantine systems, and, in particular, to applying 
to rational functions as produced by the MacMahon lifting. 

Let r E IK(xi,..., xa) be a rational function whose denominator factors into binomials. Let elimi¬ 
nate the last k of the variables x;. Due to linearity of Q^, we can, without loss of generality, restrict our 
attention to rational functions of the form 

X* 

r(x) =- 

(1 - aix'^i) ■... ■ (1 - anX'^") 

where h,ai,..., e \ jO], cti,..., q:„ e IK and the vectors b,ai,...,a„ are forward with respect to the 
given ordering of the variables x;. Now, replace r with another rational function s e Q(zi,. . .,z„,xi,. . .,x,^) 
defined as 

x^ 

s{z,x) =-. 

(1 - Zix“l) •... ■ (1 - ZnX'^") 

This is the rational function constructed by the MacMahon lifting from the system Ax ^ -b, where A is 
the matrix with the <21 as columns. It corresponds to the simplicial cone with generator matrix () and 
apex {_*j,). Therefore the methods discussed in this paper can be applied to compute a rational function 
expression for ( 5 ). 

To obtain 0^(r) from L2^(s), all that is necessary is to substitute z; = a, for all i in 0^(5). This sub¬ 
stitution is possible by construction, as long as the rational function expression obtained for (s) is 
brought into normal form first. However, as we will see below, the rational function expressions we 
usually obtain are sums of rational functions and z = o: may well be a pole of individual summands in 
this expression. However, there are methods to do this substitution efficiently, and obtain 0^(r) from 
L2j (s), without bringing (s) into normal form first. 
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There are a number of references that discuss such methods in depth, so we will not concern our¬ 
selves further with these issues in this paper and refer the interested reader to the relevant literature: 
Substitutions into sums of rational functions of the above form are an essential part of Barvinok’s algo¬ 
rithm for counting lattice points in polytopes [16,17, 13, 31]. In [15], Barvinok and Woods extend these 
techniques and develop a set of algorithms for performing a number of operations on rational function 
expression in the above form. In particular, they give an algorithm that computes the Hadamard prod¬ 
uct of two rational functions that runs in polynomial time if the number of factors in the denominators 
is fixed [15, Lemma 3.4]. This immediately gives an algorithm for computing L2>(r) for any r of the 
form given above. This algorithm runs in polynomial time, if d and n are fixed. Computing Hadamard 
products with this method has been implemented, e.g., in the barvinok package [62, 61]. In [65], Xin 
describes another algorithm that reduces the computation of (r) to solving linear Diophantine sys¬ 
tems. 


2.6. Xin’s Partial Fraction Decomposition. This subsection makes use of material introduced in Sec¬ 
tions, including symbolic cones, prim, and the Lawrence-Varchenko decomposition. Readers unfamiliar 
with these notions should skip Section 2.6 on a first reading and return after they have finished Section 3. 

In [64, 63, 65], Xin proposes several partition analysis algorithms based on iterated partial fraction 
decomposition (PFD) in the multivariate rational function field. While [63] focuses on PFD of rational 
functions with two factors in the denominator, [65] applies PFD directly to rational functions with n 
factors in the denominator. 

As it turns out, partial fraction decomposition has a clear interpretation in the language of polyhe¬ 
dral geometry: Given a simplicial cone C c IR"^, PFD writes C as an inclusion-exclusions of simplicial 
cones Ci,...,Cd, such that, for each C/, all but one of the generators of C* lie in the Xd - 0 coordinate 
hyperplane. This decomposition is closely related to a Brion/Lawrence-Varchenko decomposition of a 
section of C. 

For example, consider the simplicial cone C c IR^ with apex 0 and generators (1,0,1), (0,1,1) and 
(2,1,2). The corresponding rational function is 

1 

(1- xz){l-yz){l-x^yz'^)' 

Computing the PFD, we find 

1 xy~^ x-i- x^ -i- x^z-\- x^yz 

O f~' “ -— -“I” -• 

(1 - xy“i)(l - x^y^Kl - yz) (1 - xy“i)(l - y)(l - xz) (1 - x^y^Kl - y)(l - x^yz^]' 


On the level of symbolic cones, this corresponds precisely to the decomposition 


( 10 ) 


- K"” ((-:■ V1 - K“'((-;■ nill ((V 


Geometrically, this decomposition can be obtained as follows. The intersection of C with the hyper¬ 
plane {(jc, y, z) I z = 2} is the triangle A with vertices (2,0,2), (0,2,2) and (2,1,2). Theheight z = 2 (slowest 
possible such that the triangle has all integer vertices. A has Lawrence-Varchenko decomposition 

=K” ((-;■-■) fi))i - K” ((-:■■) 


Each of these three 2-dimensional cones, we now translate to the origin and take the Minkowski sum 
with the ray through its former apex. This means simply replacing vf); cj) with5fg'’^’'’^’°*((r'i, 

V 2 , prim(( 7 ))). Thus, we get precisely the decomposition given in (10). We have now obtained a (signed) 
decomposition [C] = ci [Ci] -i-... -i- Cd [Cd] such that each Q has only one generator with z- XdfO- Such 
Ci have the desirable property that the intersections Ci n (x | S 0} are much easier to describe. On 
the rational function level, Xin [65] makes use of this fact by using the PFD of a rational function p - 
pi-i- ...-i- Pd as a starting point and then giving a recursive procedure in the spirit of Euclidean algorithm 
formulas for computing the 0^(p;). We will come back to Xin’s algorithm in Section 5. By contrast, the 
algorithm we describe in Section 3 makes direct use of the Lawrence-Varchenko decomposition: We 
give explicit formulas on the level of symbolic cones that provide a Lawrence-Varchenko decomposition 
of the polyhedron C n {x | x^; S 0} itself. 
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3. Polyhedral Omega- the New Algorithm 

In this section, we introduce our new algorithm, Polyhedral Omega, which is a synthesis of the parti¬ 
tion analysis approach to solving linear Diophantine systems and methods from polyhedral geometry. 

The Polyhedral Omega algorithm solves the rdfLDS problem, as defined in the introduction. Given 
a matrbc A e and a vector b e Z"*, it computes a rational function expression p in d variables 

zi,...,Zrf whose multivariate Laurent series expansion represents^ the set of all non-negative integer 
solutions X £ Z^Q of Ax b. As motivated in Section 2, this rational function expression corresponds to 
a representation of the set {x e Z^q | Ax ^ b} as a signed linear combination of simplicial cones. Such 
a representation in terms of what we call symbolic cones, defined in Section 3.2 below, is central to the 
algorithm and can serve as a useful output of the algorithm in its own right. 

Without loss of generality, we will deal only with systems Ax ^ b of inequalities in this paper. How¬ 
ever, it is straightforward to extend the algorithm presented here to handle mixed systems containing 
both inequalities and equations directly. In particular, we have added this feature to our implementa¬ 
tion of Polyhedral Omega to which refer the interested reader for details [23]. 

The Polyhedral Omega algorithm consists of the following steps: 

(1) MacMahon Lifting. Given a matrix and a vector h e Z"*, we use the MacMahon lifting 

to compute a unimodular simplicial cone C c such that 

|x e J Ax = Oj(O^(■ ■■(0^(0))) = (C), 

I I > ^ ‘ V 

m times 

i.e., the desired set of solutions can be obtained by applying the Omega operator which elimi¬ 
nates the last coordinate m times to the cone C. 

(2) Iterative Elimination of the Last Coordinate using Symbolic Cones. This is the main step of 
the algorithm, in which we compute 12^ (C) iteratively, by eliminating one variable at a time. 
Each application of corresponds to intersecting a polyhedron P c R” with the half-space 

{xeK’^l Xn> O} of points with non-negative last coordinate and then projecting-away 
that last coordinate. The crucial property of this procedure is that, throughout, polyhedra are 
represented as linear combinations of simplicial symbolic cones C which are stored as triples 
(V, q, o) of a generator matrix V, an apex q and an openness-vector o. In the end, we obtain a 
representation 

0^(0 = X a; Q 

in terms of symbolic cones, much more compact than a representation in terms of rational 
functions. 

(3) Conversion to Rational Functions. The representation of the solution in terms of symbolic cones 
is then converted to a rational function representation. For this, two different methods can be 
used that each have their own advantages and that can also be used in tandem. 

(a) One option is to explicitly enumerate the lattice points in the fundamental parallelepiped 
for each of the C,- and then apply equation (2) to obtain a rational function. 

(b) Another option is to use Barvinok decompositions to represent each Q as a short signed 
sum of unimodular cones, which can be immediately converted to rational functions. 

(4) Rational Functions in Normal Form. Each of the alternatives in the previous step gives a ra¬ 
tional function expression for the solution. However, this rational function expression will not 
in general be in normal form. Standard algebraic procedures may be used to bring the ratio¬ 
nal functions on a common denominator and to simplify the resulting expression. Note that 
conversion to normal form may increase the size of the output exponentially. 

Depending on the application, it may be a good idea to work with the intermediate representations 
computed over the course of the algorithm and skip the later steps. As already mentioned, the represen¬ 
tation in terms of symbolic cones is in most cases significantly more concise than the rational function 
representations. Especially if the output of the algorithm is to be examined by hand, we recommend 

^The coefficient of := ■ ...■ z^^ in the Laurent expansion of p is 1 if x is a solution of Ax » fi, x » 0 and it is 0 otherwise. 
In the following, we will sometimes abuse terminology and refer to the coordinates of solution vectors and the corresponding 
variables of the generating function interchangeably. 
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taking a look at the symbolic cone representation first. Among the rational function representations, 
Barvinok decompositions can be exponentially shorter than expressions obtained from enumerating 
fundamental parallelepipeds. On the other hand, the enumeration of the fundamental parallelepipeds 
can be described in terms of explicit formulas and expressions with fewer distinct factors in the denom¬ 
inators of the rational functions. Computing the normal form of a rational function can lead to useful 
cancellations, but is computationally expensive and can also increase the size of the output exponen¬ 
tially. 

We now go through each of the steps of the algorithm in detail. 


3.1. MacMahon Lifting. Just as MacMahon did, we begin by using the MacMahon lifting, which we 
already met in Section 2, to transform the representation of the problem. The MacMahon lifting is 
implemented by the function macmahon defined in Algorithm 1. 

Given A e and hE Z™, we construct a matrix V e and a vector q e by 






Geometrically, we view this construction as defining a cone C = ^(V; q) with apex q and the columns 
of V as generators. As already mentioned, the crucial property of C is that 

n^iC) :=Os(ns(---(Os(C)))) = {zeK^o I Az^bl, 

m times 

where for any polyhedron P e IR”, 

and :;r : K” ^ is the projection that forgets the last coordinate. Moreover, this construction of C 
is such that for any polyhedron P obtained throughout the process of applying O^, the projection n 
maps P bijectively onto its image and moreover preserves the lattice structure of P. More precisely, let 
X denote the affine hull of {x | Ax ^ b,x^0] and let (X) denote the result of projecting-away 

the last j coordinates. Then n : ^ induces a bijection between XJ~^ and XK ft follows 

that n also induces bijections between P and n{P) for any P c X^~^. This property of the MacMahon 
lifting will be crucial for our construction below. In fact, we have the even stronger property that n 
gives a bijection between n X^~^ and n XP This latter property is essential when working 

with rational functions, but it is irrelevant for our purposes, as we will be working with symbolic cones 
instead. 


Algorithm 1: MacMahon Lifting 
Input : A matrix A e Z'"’'^ and a vector b e Z"*. 

Output: A symbolic cone C = {V, q, o) with V e ji<i+m)xd^ ^ ^ j^d+m ^ ^ s^^-h that 

D.'^[C)^{zeU^^\Az>b}. 
def macmahon(A, h): 



o = 0e{0,1}^ 
return (V, q, o) 


3.2. Iterative Elimination of the Last Coordinate using Symbolic Cones. A symbolic cone C is simply 
a triple [V, q, o) where 1/ is a matrix of generators, q is the apex of the cone and o is a vector indicating 
which faces of C are open. This is the form in which we are going to store cones throughout the algo¬ 
rithm. The triple (V, q, o) represents the set q), just as defined in Section 2. A linear combination 

of symbolic cones is a formal sum X ctiQ. i-e., a list of pairs (o:,, Q) where each Ci is a symbolic cone 
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given as a triple. The meaning of these linear combinations is best understood in terms of indicator 
functions. For any set S c we define the indicator function [S] : ^ R by 


[S](x) = | J 


if X e S, 
if X ^ S. 


Indicator functions form a vector space. In particular, addition of indicator functions is defined point- 
wise, as is multiplication with a scalar. For a symbolic cone C = (V, q, o) we will write [C] to denote the 
indicator function q)] of the associated simplicial cone. Thus, the formal sumX (tjQ represents 

the indicator function X oci [Q]. Later on we will also talk about equivalence classes of such indicator 
functions, but it is important to stress that by default [S] stands for the actual indicator function and 
not an equivalence class thereof. 

In the second step of the algorithm, we compute O J (C) by iteratively applying a function el i m LastCoord, 
summarized in Algorithm 3, that takes a linear combination Y.ociCi of symbolic cones and returns a 
representation of (X oLi [C/]) as a linear combination of symbolic cones. 

The action of the Omega operator on indicator functions is induced by its action on sets. There is, 
however, an important subtlety. The intersection with the half-space is straightforward to define for 
indicator functions: For any indicator function /, 




if X e H^, 
ifx^H+. 


To define the projection, we make use of the property of the MacMahon lifting that for any set S appear¬ 
ing throughout the algorithm and any (xi,. ..,x„-i) e Ji{S) there exists a unique Xn such that (xi, ...,x„) e 
S and the same holds for indicator functions. This allows us to define n{f) for an indicator function / 
by 




/(Xi,...,x„_i,x„) if 3x„ : /(Xi,...,x„_i,x„) ^ 0 
0 otherwise. 


Since, if it exists, the x„ in the above definition is uniquely determined, it follows that n{f) is well- 
defined. Actually, all indicator functions / to which we apply n moreover satisfy f{x\,..., x„_i, x„) = 1, 
however, we will not make use of this fact. With these preparations we can now define 


ns(/):=;r(H+n/). 


Since the Omega operator is linear, we have the important consequence that 

ns (ya/[Q]) = ya;n^([Q]). 

Therefore it suffices to define elimLastCoord on symbolic cones Q - we do not need to handle general 
polyhedra or indicator functions directly. 

To compute (C) for a simplicial symbolic cone C, we decompose n C into a signed sum of 
cones. The most well-known decomposition of this type is the Brion decomposition which states, 
roughly, that “modulo lines” a polyhedron is equal to the sum of the tangent cones at its vertices. Let us 
make this statement precise. For a given polyhedron P and a point v, the tangent cone of P at u is 

tconelP v) = {v+ u \ v + due P for all sufficiently small d > 0 }. 

Similarly, the feasible cone of P at f is 

fconelP v + due P for all sufficiently small d > 0 }, 

i.e., tconelP v) - v + fconefP v). We call the tangent cones of P at its vertices v the vertex cones of P. For 
each vertex cone C/ = tconefP vi) of P let pc, denote a rational function expression for the generating 
function (Pc, of lattice points in C*. For simplicial cones, pc, can be obtained via (9) as we have seen in 
Section 2. For non-simplicial cones, such expressions can be obtained via triangulation, for example, 
as in Figure 6. Fortunately, it turns out that we will only have to deal with the simplicial case in our 
algorithm - we never have to triangulate. Given the above notation, Brion’s theorem states the following. 


Theorem 3.1 (Brion [24]). For any polyhedron P with vertices f i,..., ujv, 

N 


®p(x) — V. PtconefPi^.l(x). 

i=l 
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Geometrically, i.e., on the level of linear combinations of indicator functions, this identity can be 
read as 


N 

[P] = ^ [tcone(_fi Vi)] modulo lines . 

j=i 

The phrase “modulo lines” means that this identity holds true up to a finite sum of indicator functions 
of polyhedra which contain lines. A line is a set of the form a+Rb for some vectors a, b. See [13] for 
details, in particular Theorem 6.4 therein. 

For our purposes, working modulo lines is not ideal, since the Omega operator is not well-defined 
modulo lines, as we already discussed in Section 2.4. Another example that shows this, even in the 
simple case of one-dimensional intervals^, is this; 


n^((-oo, 1]) = [0,1] ^ -(1,00) = -ns((i,oo)) 


even though (-oo, 1] = -(l,oo) modulo lines. We therefore have to work with a decomposition that 
holds exactly, not just modulo lines. 

One option is the Lawrence-Varchenko decomposition, which makes use of the notion of forward 
cones [45]. In general, “forward” can be defined with respect to an arbitrary half-open half-space. For 
our purposes, the following more particular definition will suffice, however: A vector v is forward if its 
first non-zero entry is positive and a symbolic cone is forward if all of its generators Vi are forward. Note 
that the cone obtained from the MacMahon lifting is forward by construction. 

For any non-zero vector v, we have that v is forward if and only if - y is not forward. Therefore, we 
can turn a simplicial cone into a forward simplicial cone by reversing the direction of all “backward” 
generators - as long as we adjust the sign and which faces of the new cone are open. The resulting 
forward cone will be equivalent modulo lines to the original cone. Here we make use of the fact that for 
any vector vf^O, the one-dimensional ray'^°(y) satisfies ['^^{v)] = v)] modulo lines. 

In general, let C = {V, q, o) denote a symbolic cone, with V e q and o e [0,1}*^. Let bwd(C) 
denote the number of columns of V that are backward, i.e., not forward, and define sgn(C) (- 
Let V denote the matrix obtained from V by multiplying all backward columns by -1. Let o' denote the 
vector obtained from o by letting o'j- I- oj if the j-th column has been reversed and oi = oj otherwise. 
Define flip(C) := {V',q,o'), as summarized in Algorithm 2. Then we have that flip(C) is forward and 
moreover 

[flip(C)] = sgn(C)[C] modulo lines. 

Applying this operation to Brion’s theorem, we can obtain the following identity of indicator functions, 
which is essentially the theorem of Lawrence-Varchenko, even though their theorem was originally 
stated in terms of generating functions. 


Algorithm 2: Flip 

Input : A symbolic cone C = (V, q, o) with V e qeQ'^ and o e [0,1}*^. 

Output: A pair (s, C') of a symbolic cone C' and a sign s such that [C] = s[C'] modulo lines and 
s = sgn(C). 
defflip(C): 

/(() = the i-th column of V is backward 
o{i) - 1 if/(/) else -1 

V' - the matrix obtained from V by multiplying the i-th column by a{i) 
o' - the vector defined by o'. = xor(oj,/(i)) for i- 1,..., fc 
C'^{V',q,o') 

return (s, C') 


3 


Here, we simplify notation and write [a,b] instead of [[a,,^]] for the indicator function ofthe interval [a,b]. 
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Theorem 3.2 (Lawrence-Varchenko [45, 60]). For any polyhedron P with vertices Vi,..., 

N 

[P] = ^ sgn(tcone(P, i^/))[flip(tcone(P vt))] 
i=l 

The crucial point here is that this identity holds exactly, without taking equivalence classes modulo 
lines. The advantage of this exact identity over Brion’s result is that we can use it safely in conjunction 
with the Omega operator. Since we are not taking equivalence classes, the problem of well-definedness 
does not even arise. Working with rational function representations necessitates working with equiva¬ 
lence classes, a complication we can avoid by working on the level of symbolic cones instead. 

The key observation behind the Polyhedral Omega algorithm is now that for the simplicial cones C 
appearing through the course of the algorithm, we can give explicit formulas for a Lawerence-Varchenko 
decomposition of the polyhedron C n into simplicial cones. This elimination of the last variable is 
performed by the function el imLastCoord given in Algorithm 3. To see the correctness of this algorithm, 
we proceed just as we did above; We first derive explicit formulas that determine a Brion decomposition 
and then we convert this Brion decomposition into a Lawrence-Varchenko decomposition by flipping 
all cones forward. Since the resulting decomposition is an exact identity of indicator functions, the 
Omega operator can safely be applied to the resulting decomposition in the next iteration. 


Algorithm 3: Eliminate Last Coordinate 


Input : A symbolic cone C = 


(V, q,o) with V e q e < 


and o E { 0 , 1 }*^, such that the 


projection n that forgets the last coordinate induces a bijection between aff(C) and 
Output: A linear combination Z / o: i Q of symbolic cones, represented as a dictionary mapping 
cones Ci = (Vi,qi, o;) to multiplicities a;, such that 0^([C]) = n/IC;]. 

def elimLastCoord(C): 

Vj - the 7 -th column of V 

r = {; e ik] I vj,n > 0} 

J-^ ljelk] I Vj,n<0} 

o'( 7 ) = (o'( 7 );),=i,...,i; where o'{j)i = o, if i 7 else 0 


i-v 


T-i = 


],n 


Vl.n 


^j,n 


-1 


^ j+l,n 
-Vj,n 


^k,n 


V 

GJ = prim(:/T(sg(q„) ■ V-tJ)] 

B+ = for jef^] 

B~ = [(G-'’,7r(w7),o'(7)) for /e/"] u [{pnmin{V)),Ji{q),o)] 

B - B~ if qn^ 0 else B^ 

D - the dictionary that maps flip(C') to sgn(C') for all C' e B 

_ return D 


To derive these formulas, we have to distinguish three cases: Whether the apex <7 of C lies above the 
hyperplane = {x | = 0 }, on the hyperplane or below the hyperplane. 

Note that our original matrix A has m rows and d columns. After the MacMahon lifting, we are then 
dealing with a d-dimensional cone C in Through successive eliminations of the last coordinate, 

we will eventually obtain a linear combination of d-dimensional cones in K'^. Throughout this process 
we will denote the index of the last coordinate with n-d+m,d+m-\,...,d and we denote the number 
of generators by fc = d. 
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(A) Intersecting a cone with the half-space with 
apex above the hyperplane 



(b) Brion decomposition of the intersection. 


Figure 11. Intersecting a cone with the half-space when the apex is above the 
hyperplane 


Case A: Apex above the hyperplane. Let C be a simplicial cone in K" with apex q, openness o and gen¬ 
erators vi,...,Vk such that t/n > 0. This situation is illustrated in Figure 11. Since q is contained in the 
half-space the polyhedron C n has q as one of its vertices. The other vertices of C n are 
in one-to-one correspondence with the intersections wj of the extreme rays Rj of C with the hyper¬ 
plane Hn- Since qn > 0, the extreme ray Rj generated by Vj intersects Hn if and only if vj^n < 0, i.e., the 
generator vj points “down” with respect to the last coordinate. Let 

r = {;■ e [fc] I vp„ < 0} 


denote the set of all indices j such that the corresponding generator vj points down. The intersection 
Wj - HnH Rj of the intersection hyperplane with such a “downward” ray Rj, j ^ J~ can be computed by 

Wj - q - Vj. 

Vj,n 


Applying Brion’s theorem, we find that, modulo lines, C n can be written as a sum of vertex cones 
at the vertices q and Wj for j e }~. The vertex cone at q is identical to C, in particular it has the same 
generators. The generators ,..., of the vertex cone Cj tcone{Wj,C n fall into three classes: 

(1) One of the generators of Cj, let us call it g^, is just the reverse of the generator Vj of C, i.e., 

■' 7 J 



In other words, gj points from the apex Wj of Cj to the vertex q of Cn H^. 

(2) Foreach/E /“within 7 , there is a generator g^^ of Cj that points to another vertex w,-of CnHjj. 
These generators gj are given by 

/ ^ j.n j ^ Vj^nVj-Vj^nVi- 

qn 

(3) Finally, for each i e [A:] \ /“, there is a generator gj that lies in the intersection of Hn with the 
two-dimensional face of C that is generated by Vj and Vj. Note that Vj points down while Vj 
points up or is horizontal with respect to the last coordinate. We can therefore find a suitable 
gj by taking a linear combination of Vj and Vj with non-negative coefficients, such that the last 
coordinate of gj is zero. Thus, we obtain 

gj = Vi.nVj-Vj^nVj. 

If Vj is horizontal, i.e., Vj^n - 0 . then gj is just a multiple of Vj. 
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Note that even though, intuitively, the generators of the second and third kind arise from a different 
construction, the formulas defining gj in both cases are identical; gj is a linear combination of vj and 
Vi such that the last coordinate of gj is zero. 

It is important to note that the cones tconeCy, C n for v a vertex of C n are all simplicial. If 
V - q, then this is true by assumption. If h = wj for some j, then we can observe this by noting that 
the matrix G^, which has the generators gj of Cj as columns, arises from the matrix V, which has the 
generators vt of C as columns, via the transformation -V -T where T is a.kx k matrix of the form 


( 


\ 




V\,n 


^j,n 

Vj-l.n 


-1 


^ i+l,n 
- V],n 


^k,n 


V ^j,nJ 

which is of full rank as vj^n < 0 . 

We have now determined apex and generators for each of the cones Cj. Moreover, we observe that 
the half-space H'^ is closed. Therefore the facet of Cj corresponding to generator gj is closed, while all 
other facets inherit their openness from C. The openness vector o'{j) = (o' ( 7 ) i,..., o' (j) fc) of Cj is thus 
given by o'( 7 ) 7 = 0 and o'[j)i = oi for iV 7 . 

We thus have found an explicit Brion decomposition 

[Cn H+] = [C] -h Y. IQ] modulo lines, 
jer 


To turn this into an explicit Lawrence-Varchenko decomposition of Q>(C), all we have to do is to a) 
flip all cones forward, by applying the map flip defined above, and b) project away the last coordinate, 
by applying the projection n : (jci,..., x„) ^ (xi,..., x„_i). Both of these operations can be carried out 
easily on the level of symbolic cones. Here it is important to recall two important properties: On the 
one hand, flip may introduce a sign in front of each cone and may change the openness of a cone by 
toggling the openness of each flipped generator. On the other hand, the affine hull of aff(C) of C is such 
that n^laff(C) is a bijection onto its image, so that in particular Ti\cj is one-to-one for each cone Cj in the 
decomposition. With these observations we obtain a decomposition 

[Os(C)] = [;r(CnHQ] = [;r(C)]+ Y sgnfCyjlyrCflipfCj))] 


which holds exactly, i.e., not just modulo lines. Note that C does not need to be flipped since we know 
that C is forward, by construction. 


Case B: Apex below the hyperplane. The case where the apex q of cone C lies strictly below the hyper¬ 
plane Hn, i.e., qn < 0, is similar to the previous one. The difference is that here, q is not a vertex of 
C n as it is cut off by the intersection, as illustrated in Figure 12. The vertices wj of Cn are thus 
in one-to-one correspondence with generators of vj of C that point “up”, i.e., that satisfy vj^n > 0. Let 

/+ = {7 E [k] I Vj,n > 0} 

and let wj n Rj denote the intersection of with the ray Rj - q + IR^o fj generated by vj for each 

7 E /^. Then w; can be computed via 

qn 

Wj - q - Vj. 

Vj,n 

Let Cj = tcone(W 7 ,C n H^) denote the vertex cone at wj. Again, we can think of the generators 
gj,...,gj of Cj as belonging to three different types. 

(1) One generator, which we call gj, is just the corresponding generator of C, i.e. 
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(A) Intersecting a cone with the half-space (B) Brion decomposition of the intersection, 

when the apex is below the hyperplane . 

Figure 12. Intersecting a cone with the half-space H'^ when the apex is helow the 
hyperplane 


(2) For every ie J^,i^ j, there is one generator gl which points from wj to the vertex WiofCn , 
i.e. 


si = 


^j,n ^i,n 


{Wi - Wj] = Vj^n Vi - Vi^n I’j- 


Note that the sign is opposite from what we had in the previous case to ensure that the factor 
_ Vj,nV, n positive, i.e., gj is a positive multiple of wi - wj. 

(3) For every i e [k] \ , there is one generator gl which lies in the intersection of with the face 

of C generated by Vj and y,-. Again, we take gj to be a non-negative linear combination of Vi 
and Vj, namely 

gj = Vj,nVi - Vi,nVj- 

If Vi is horizontal, i.e., = 0, then gj is just a multiple of Vi . Again, the terms in this expression 

have opposite signs compared to the previous case, in order to guarantee that the coefficients 
Vj fi and - Vi ji are both non-negative. 

Just as in the previous case, the expressions for generators of the second and third kind are identical. 
Also, we can set up a transformation matrix T' to express the generator matrix G-' of Cj as a transforma¬ 
tion of the generator matrix V of C, via = V ■ T'. In this case, T' is given by 


Vj,n 


I 


r' = 


-Vl,n 


Vj,n 

~Vj-l,n 1 


^j+l,n 


^j,n 


^k,n 


V 


Vj,n 


First of all, we note that this has again full rank since Vj^„ > 0, which guarantees that the cones Cj are 
simplicial. But more importantly, T' - - T, which means that both cases A and B can be treated with 
one common transformation 

& = sgiqn]-V-T. 

where 

sg(x) = 1 if X S 0 else - 1. 

The definition sg(0) = 1 will come in handy in case C, as we will see below. 

We now have computed generators and apices for the vertex cones Cj. Again, the half-space is 
closed, whence the facet of Cj corresponding to generator gj is closed, while all other facets inherit 









POLYHEDRAL OMEGA: SOLVING LINEAR DIOPHANTINE SYSTEMS 


27 



V3 

(A) Intersecting a cone with the half-space 
when the apex is on the hyperplane 



(b) a decomposition of the intersection is a given by 
a deformation of the Brion decomposition in case A. 


Figure 13. Intersecting a cone with the half-space when the apex is on the hy¬ 
perplane Hn- 


their openness from C. Therefore we now have an explicit Brion decomposition 

[C n H'^] = ^ [Cy] modulo lines. 
j^r 

We proceed just as in the previous case and apply flip and n to obtain a Lawrence-Varchenko decompo¬ 
sition 

[QslO] = [:/r(CnH^)] = ^ sgn(Cy) [7t(flip(Cy))] 

that holds exactly, not just modulo lines. 

Case C: Apex on the hyperplane. At first glance, the case where the apex q of cone C lies on the hyper¬ 
plane Hn, i.e., qn = 0, appears special. The reason is that in this case C n can be a non-simplicial 
cone, as shown in Figure 13(a), and that the only vertex cone of C n f/jj is C n itself Thus a simple 
Brion decomposition would take us out of the regime of formal sums of simplicial cones. Of course we 
could solve this issue by triangulating C n But triangulation algorithms are a complex subject in 
their own right [34,46, 53], at odds with the simple approach based on explicit decomposition formulas 
we have applied thus far. 

Fortunately, it turns out that the same formula we used to handle case A can also be applied in case 
C to obtain a decomposition into simplicial cones. It just so happens that all simplicial cones that 
appear have the same apex q. This phenomenon can be explained with a deformation argument; For 
a small e > 0, the apex of ce„ -i- C lies above the hyperplane Hn so that case A applies and we obtain a 
decomposition 

[(ee„ -t C) n H^] = [een -i- C] -i- ^ [Cy (e)]. 

HI 

The generators of the cones (ee„ -i- C) and Cj (e) are independent of e. The apices iee„ + q] and Wj (c) 
of the cones (ce„ -i- C) and Cy(e), on the other hand, do depend on e. However, this dependence is 
continuous so that (ee„ + q] ^ q and wj (e) —► as c ^ 0. In the limit, we thus obtain the desired 
decomposition of Cn Hn, as illustrated in Figure 13(b). That this works in general is the content of Liu’s 
theorem. 

Theorem 3.3 (Liu [48]). Let Ae IR'”’'”, h ; (-1,1) ^ IR"* a continuous function andP{t) - {x \ Ax ^ b{t)]. 
Suppose that for each 0 ^ te (-1,1) the polyhedron P{t) is non-empty and has precisely £ vertices Wty,..., 

. Suppose further that for each j - there exists a fixed cone Cj such that 

fcone(w t,j,P{t)) - Cj forallO^ te (-1,1), 

i.e., the apex of each tangent cone tcone{wt,j,P(t)) may change depending on t, but the feasible cone 
fcone{wt,j,P{t)] does not. Then, for each j, the vertex Wtj converges to some vertex ofP{0) as f —► 0. Let 
V be a vertex o/P(0) and let Jy be the set of all j such that Wt,j converges to v. Then 

[tconelf,P(0))] = [v-tCj] modulo lines. 

HJv 
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To apply this theorem, we first observe that there exists a matrix A and a continuous function h 
defined on the half-open interval [0,1) such that 

(ee„-I-C)n -{x\Ax^ b(e]}. 

To extend b continuously to (-1,1) we simply let h(-e) b[e). By construction, {ecn -i- C) n is non¬ 

empty for e > 0. Moreover (ee„ -i- C) n H'^ has the same combinatorial type for each e > 0. In particular, 
the number of vertices of (ee„ -i- C) n and the feasible cones at these vertices do not change for e > 0. 
For e = 0 we obtain the polyhedron C n H'^ that we are interested in. C n has just one vertex, namely 
q, to which all the vertices of (ee„ -i- C) n converge. Therefore, the conclusion of Theorem 3.3 gives 
us the desired decomposition 

[Cn H+] = [C] -h Y. tC'fl modulo lines 

Hr 

where each Cj is given by precisely the same formulas as in case A. In particular, the formula for Wj 
gives q as the apex of each of the Cj, as expected. We can view this decomposition as giving, implicitly, 
a triangulation of the dual cone of C n lust as in the previous two cases, we obtain an exact de¬ 
composition of C n (not just modulo lines) by flipping all cones on the right-hand side forward. By 
projection we obtain 

[Os(C)] = [;r(CnH+)] = [;r(C)]+ Y sgn(C;)[;r(flip(Cj))]. 

Hr 

Note that our convention sg(0) = 1 ensures that = sg(^ 7 „) -V -T holds for all three cases. 

Summary and Optimizations. We have now seen that elimLastCoord correctly applies to a linear 
combination of symbolic cones of the form produced by macmahon. The function elim, defined in 
Algorithm 4, iterates this elimination of the last variable m times. As the output of elim we thus obtain a 
list of simplicial symbolic cones Ci,..., Cn and multiplicities a i,..., such that 

f I 1 ^ 

( 11 ) = 

i=i 

where each symbolic cone C; is given in terms of a triple (V, o, q) where V is an integer matrix of gener¬ 
ators, o is a vector indicating which faces of the cone are open and q isa rational vector giving the apex 
of the cone. Note that each C; is a d-dimensional cone in The above identity of indicator functions 
holds exactly, not just modulo lines. 


Algorithm 4: Eliminate Coordinates 

Input : A symbolic cone C-[V,q,o) with V e ^ ^ Qd+m ^ ^ 2 }^ such that 

projection that forgets the last m coordinates induces a bijection between aff(C) and 
Output : A linear combination T-i ai Ci of symbolic cones, represented as a dictionary mapping 
cones Ci - {Vi,qi, Oi) to multiplicities a/, such that 02*([C]) = Y.i tt/ICfl. 
defelim(C): 

for i-\,...,m 

|_ C ^ map(elimLastCoord,C) 

return C 


The elimination algorithm is easy to implement, since all cones are stored in terms of integer matri¬ 
ces and all transformations are given explicitly. However, it is important to highlight two implementa¬ 
tion aspects for performance reasons. 

First, it is important to ensure that the generators for a given cone C,- are stored in primitive form, 
as defined in Section 2. Algorithm 5 shows how to compute prirnff) for integer vectors. Replacing a 
column V of V with primff) does not change the cone C; the generator matrix V defines, but it may 
significantly reduce the encoding size of the matrix V. Therefore, after each round of elimination, all 
columns v of V should be replaced by primfy), respectively, since the transformations applied in each 
round may introduce non-primitive columns. 
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Algorithm 5: Make a vector primitive 
Input : Anon-zero integer vector v e Z”. 

Output: The shortest integer vector that is a positive multiple of v. 
def prim(i'): 

g = gcd(Hi,...,H„) 

return 


Second, one should note that after each round of elimination, a given cone C/ may appear multi¬ 
ple times in the sum (11), with multiplicities of opposite signs. It is therefore crucial to collect terms 
after each elimination and sum up the multiplicities of all instances of a given cone Q, as a signifi¬ 
cant amount of cancellation may occur in some cases. For this purpose, the generators of the cones 
C should he stored in (lexicographically) sorted order, since when generators are primitive and sorted, 
every rational simplicial cone has a unique representation as a symbolic cone. Also, it may be expedient 
to store the sum (11) as a dictionary mapping symbolic cones C; to multiplicities ai throughout the 
computation, see the definition of map (Algorithm 6). 


Algorithm 6: Map Linear Combination of Symbolic Cones 

Input : A function / which maps a symbolic cones to a linear combination of cones, and a linear 
combination X; oqCi of symbolic cones, represented as a dictionary D mapping 
symbolic cones C; to multiplicities 

Output: A dictionary representing the linear combination «r/(Q) = Hij oLi Pij Cij of symbolic 
cones where /(Cj) = Coefficients are collected by the Ctj, i.e., if 

C := Cjj = ... = then the dictionary contains only the key-value-pair 
(C —>■ 

def map(D,/): 

L = an empty dictionary with default value 0 
for {C ^ a) e D 
D’ = /(C) 
for (C' p)eD' 

L L[C']^L[C'] + a-p 

return L 


3.3. Conversion to Rational Functions. The decomposition (11) into symbolic cones that we have so 
far obtained has many advantages over a rational function expression for <Dp. In particular, (11) is very 
compact, highly structured and has a clear geometric meaning, which makes it suitable for further pro¬ 
cessing by human and machine alike. Therefore, when applying our algorithm, we recommend inves¬ 
tigating whether the representation (11) can be used instead of a rational function expression. 

In case a rational function expression for Op is indeed required, we present two different methods for 
computing it in this section. The first method explicitly enumerates the lattice points in the fundamen¬ 
tal parallelepiped of a cone, using an explicit formula based on the Smith normal form of the matrix of 
generators. The second method uses Barvinok’s recursive decomposition to express a given cone as a 
short signed sum of unimodular cones. The first approach has the advantage of giving the rational func¬ 
tion explicidy in terms of a simple formula, rather than appealing to a complex recursive algorithm. The 
second approach has the advantage that it produces much shorter formulas. In fact, if the dimension 
is fixed, Barvinok’s algorithm is guaranteed to produce polynomial size expressions, whereas the num¬ 
ber of points in the fundamental parallelepipeds involved may very well be exponential. However, the 
representation (11) in terms of symbolic cones is shorter than either of the two rational function expres¬ 
sions. Moreover, it turns out that in practice, conversion to rational functions is usually the botdeneck 
of the algorithm, taking much longer than the computation of (11) on a large class of inputs, no matter 
which of the two conversion methods is used, see Section 4. Both approaches have in common that 
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Figure 14. The Smith normal form can he viewed as aligning the bases of Z" and 
a suhlattice L so that the fundamental parallelepiped of the sublattice is a rectangle 
with respect to the chosen bases. In this example, the initial basis of L is given by the 
columns of 1^ = (^2 2 )- The Smith normal form is 1/ = USW with 17 = j), IV = (J f) 

and S = (0 g). Note that U~^ = (} j) and W~^ - 

they convert (11) into a rational function identity one simplicial cone at a time. If pc, [z] is a rational 
function expression for <l)c, iz] for any i, then 

N 

Therefore, we can restrict our attention to computing rational function representations pc for a given 
simplicial symbolic cone C. 

Fundamental Parallelepiped Enumeration. As we have already seen in Section 2, 

, ^ XpeZ”nn°(Pi,...,y,i;q) - 2 *^ 

Pc(z) - - 

where the Vi are the generators of C and n" (fi,..., : ^7) is the fundamental parallelepiped of C with 

given openness o and apex q. The number of lattice points in the fundamental parallelepiped 11 is 
equal to the determinant of the matrix V of generators (provided V is square) and this determinant is 
exponential in the encoding length ofV even if the dimensions of V are fixed. Therefore the polynomial 
in the numerator of the above expression can be exponential in size, even in fixed dimension. 

Nonetheless, this representation of pc is attractive because the numerator is very simple to compute. 
Contrary to what one might think at first glance, enumerating the lattice points in the fundamental 
parallelepiped does not require solving another linear Diophantine system. Instead, this set of points 
can be generated directly, given the Smith normal form of the matrix V. This technique is well-known 
[26, 42], but the method is very instructive and the formula we use is slightly different from [26, 42] so 
we provide its derivation here. To simplify the exposition, we will assume that the apex <7 = 0 is at the 
origin and that o = 0, i.e., the cone C is closed. 

Every integer matrix V e Z'"’'” has a Smith normal form which we define as a triple of integer ma¬ 
trices U,S,W with V - USW such that the following properties hold: First, [/ e zand W e Z”’*” 
are unimodular, i.e., they are invertible over the integers or, equivalently, det(7 = detlV = +1. Sec¬ 
ond, S e z"”*” is a diagonal matrix with diagonal entries Si,...,Sjnm(m,n] such that if A: = rankV then 
Si,...,Sfc > 0 and Sfc+i,...,Smin(m,n) = 0 and, moreover, Si| 52 , S 2 IS 3 , ..., Sk-i\sjc. The matrix S is uniquely 
determined, but the transformations U and W are not. Nonetheless we will simply say that USW - V 
is “the” Smith normal form of V. Note that if V is square then det V = si • S 2 •... ■ s^. Since in our situation 
C is always a simplicial cone, V will always have full rank so that k is the number of generators of C. For 
now, we will make the additional assumption that C is full-dimensional which means that V is square 
and k - m - n. 

One way to interpret the Smith normal form is the following. Let L denote the lattice spanned by 
the columns of V. L is a sublattice of the integer lattice Z”. In this setting, the transformations U 
and W can be viewed as changing bases on both L and Z” so that with respect to the new bases the 
fundamental parallelepiped of L is just a rectangle with side-lengths given by the diagonal elements Sj. 
Figure 14 illustrates this idea with an example. Suppose we start out with the standard basis ei,...,e„ 
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of Z" and with a given basis v\,...,Vn e Z” of L. The vi are the columns of V or, in our application, 
the generators of a simplicial cone C. The matrix U now transforms the basis of Z” into a new basis 
fip.e'„ where e'. = Ujjej. Simultaneously, W transforms the basis of L into a new basis 
where f'. = ITTl y,-. S now gives the coordinates of the relative to the basis ep...,e'„. The fact 
that S is diagonal now means that f'. is simply a multiple of e'., or, in other words, that, with respect 
to the basis e^,..., e'„, the fundamental parallelepiped fllfp..., v'^) of the basis of L is simply 

a rectangle with sidelengths given by the s;. Thus, the set of lattice points in n(fp..., v'^) is extremely 
simple: it is (a linear transformation of) a Cartesian product 







u ([0,Si)z X ■■■ X [0,s„)z) 


where the intervals [0, Si)z are taken to denote the integer values in the given range. 

We now understand ffCyp..., v'^) very well, but what is the relationship between IfCi’p..., v'^) and 
the fundamental parallelepiped Vn) that we set out to enumerate? It turns out we can eas¬ 

ily transform the set of lattice points in one into the set of lattice points in the other through linear 
transformations and simple modular arithmetic, as shovm in Figure 15. Since the matrices U and W 
are unimodular, it is clear that both flli/p..., yj^) and Iflvi,..., Vn) contain the same number of lattice 
points. Moreover, for both of the two fundamental parallelepipeds If it is true that for every x e Z” there 
exists a unique point y e Z”nn such that has x-yeL. Therefore, the function /: Z” ^ Z”nn(i;i,..., Vn) 
that maps a lattice point x to the corresponding y e Z" n flli'i,..., v^] induces abijection 


between the set that we can describe and the set we want to describe. As it turns out, / has a very sim¬ 
ple explicit description: to compute f{x) we simply take coordinates with respect to the basis v\,...,Vn, 
then we take the vector of fractional parts, and then undo the coordinate transformation. More pre¬ 
cisely, for any a e IR let fract(o:) = a - [a\ denote the fractional part and for any vector x e K" let 
fract(x) = (fract(Xj))i=i . denote the vector of fractional parts. Then 

/(X) = Infract! V^Cx)). 


Putting these two observations together we obtain 

Z”nn(i^i,...,y„) = I/fract(y-i(I7([0,5i)zx...x [0,5„)z)) 

= VfractlW^S”^ ([0,5i)z x ... x [0,5„)z)) 
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Note that S ^ is the diagonal matrbc with diagonal entries ^. All other matrices are integer. For compu- 
tational efficiency it is expedient to work with integer values throughout. The common denominator of 
all values in is simply Let s'. ^ £ Z and let S' denote the diagonal matrix with entries Sj,..., s'„. 

Then we can rewrite the previous identity to 

(12) Z"nn(r'i,...,y„) = —l/(kL“^S'([0,si)z x ... x [0,s„)z) mod s„) 


where the modulus is taken componentwise, i.e., x mod a {xi mod a)!=i,...,n for a vector jc = 

The computation (12) can he performed entirely in the integers. In particular, the final division hy is 
guaranteed to produce only integer values. 

Translating (12) into the language of generating functions, we get the following explicit rational func¬ 
tion formula pc(-z) for the generating function of a symbolic cone C: 


(13) 


Pc(z) = 


—1 


mod s„) 


(l-z‘'i)-...-(l-z''«) 


The formulas (12) and (13) give a very concise and simple solution to the problem of listing lattice 
points in fundamental parallelepipeds and thus for converting cones to rational function expressions. 
The drawback of the rational function representation (13) is that the polynomial in the numerator may 
be exponential in the encoding size of C, even if the dimension is fixed. Nonetheless, computing (13) 
is much more efficient than solving a general linear Diophantine system, in the following sense; the 
Smith normal form (SNF) of V can be computed in polynomial time, even if the dimension is not fixed, 
and once the SNF of V is known, (12) provides an explicit parametrization of the lattice points in the 
fundamental parallelepiped. In particular, this gives an output-sensitive polynomial time algorithm for 
listing the points in the fundamental parallelepiped. 

As already mentioned, this method is well-known. The expression (12) is essentially the same as [42, 
Lemma 11], the main difference being that (12) avoids rational arithmetic for performance reasons. [26] 
gives a similar method based on the Fiermite normal form. 

The approach presented above can be extended to accommodate general symbolic cones C. In 
particular, it is possible to handle cones with open facets, cones with rational apices and non-full- 
dimensional cones. The fully general statement of (12) and (13) is given in the following theorem. To 
deal with open facets we introduce the notation 


a mod° b - 


a mod b 
o-b 


if a mod b^O 
if a mod b-0 


for fl, h £ Z and o £ [0,1}, andlet JC mod” b:- (x,- mod°‘ be defined componentwise for vectors 

X £ Z” and o £ [0,1}”. Similarly, [xj := ([xj),- is defined componentwise. 


Theorem 3.4. Let C be a k-dimensional simplicial symbolic cone in given by a generator matrix 
V £ an apex q and a vector o £ [0,1}^ that encodes which facets ofC are open. Let USW - V 
denote the Smith normal form ofC, so that in particular U £ Z”’'”, W £ are unimodular and S £ 
Z”’'^ is diagonal. Denote the diagonal entries ofS £ Z”’'^ hy si,..., Sfc £ Z and define s'. = ^ £ Z. Let S' 
denote the [k x n]-matrix with diagonal entries s'.. 

If there is no lattice point in the affine hull of C, then 1.'^ ^11° [V] - 0 and<t>c(z) - 0. Otherwise, letp be 
such a lattice point in the affine hull of C Abbreviate q = U~^ {q - p) e and define q""^ = [-W~^S'q\ £ 
Z*' and -W~^S'q- ^'"*£ Q*'. Let o' be defined by o'j - oj if^'^- 0 and o'. - 0 otherwise. 

Then the lattice points in the fundamental parallelepiped ofC are given by 


n'^[V;q) = { y(x) I x£ {0,...,si - 1} X ... X {0,...,Sfc-1} } where 

y(x) = —V[[W~^S'x-i- q*"*) mod"' sO -t — q 

^k 


"^Both this decision and the computation of p, if it exists, can be readily performed using the Smith normal form. 
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Figure 16. Barvinok decomposition. 1])] = ?))] - 1])] 


and the generating function ^ciz) ofC is given by the rational function expression 


i'cU) 






Theorem 3.4 gives rise to the function enumFundPar, defined in Algorithm 7, for enumerating fun¬ 
damental parallelepipeds of symbolic cones. 


Algorithm 7: Enumerate Fundamental Parallelepiped 
Input : A symbolic cone C = [V, q,o] with V e and o e {0,1}^. 

Output: The list of all lattice points in the fundamental parallelepiped of C, in no particular order, 
def en urn Fund Par(C): 

if aff(C) does not contain a lattice point 

L return [] 
else 

p = a lattice point in aff(C) 

U,S,W - SmithNormalForm( V) 
s\,...,Sfc- diagonal elements of S 
s’. = — forall i = 1,...,A: 

I Si 

S' - the (fc X u)-matrix with diagonal entries s'. 
q^ U~^[q-p) 

^frac^_p^-ly^_^int 

o' = (o^)y=i_ where o’j = oj if qj^^ = 0 else 0 
P - {(xi,.. ..Xfc) I Xj = 0,..., S; - 1 for all i = 1,..., fc} 

L = IvdlV'iS'x-i- q'°^] mod"' s^) + + Sfcr/j for xe P] 

L return L 


Barvinok decomposition. An alternative way to convert a simplicial symbolic cone C to a rational func¬ 
tion is to construct a signed decomposition of C into unimodular cones C/. A cone C; is unimodular if 
it contains just a single lattice point in its fundamental parallelepiped, which implies that the rational 
function representation (9) has just a single monomial in its numerator. 

It is crucial that we are looking for a signed decomposition of C into unimodular cones C,, i.e., for 
a way of representing C as an inclusion-exclusion of the C/. While it is always possible to write C as a 
union of unimodular cones Ct, the number of such cones Ct required may well be exponential in the 
encoding size of C, as the left-hand side of Figure 16 shows: The cone generated by (1,0) and (1, a) has 
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encoding size loga. While it can be written as a union of unimodular cones Cj generated by (1, i) and 
(1, i + 1) for i - 0,...,a- I, this positive decomposition requires a cones, which is exponential in loga. 
Signed decompositions, on the other hand, allow us to write C as a difference of just two unimodular 
cones, as shown on the right in Figure 16. If Ci is generated by (1,0) and (0,1) and C 2 is generated by 
(0,1) and (1, a) with openness-vector o = (1,0), then [C] = [Ci] - [C 2 ]. 

Barvinok’s theorem [16], which is one of the landmark achievements in the field in the last decades, 
states that this works in any dimension; Every simplicial cone C can be written as a signed sum of 
half-open unimodular cones Q such that the number of cones C* is bounded by a polynomial in the 
encoding size of C, provided that the dimension of C is fixed. That is, the number of cones will grow 
exponentially with the dimension of the matrix V of generators of C, but once the dimensions of V are 
fixed, the number of cones depends only polynomially on the encoding length of the numbers in V. 
Such a decomposition is called a Barvinok decomposition. 

The construction of Barvinok decompositions is out of scope of this article. We refer the interested 
reader to the textbooks [13, 31] for details. Briefly, the basic idea is the following. Given a simpli¬ 
cial d-dimensional cone C with generators use the LLL algorithm to And a vector w that is 

“short” with respect to a certain basis. Then decompose C into cones Ci,..., where C; is generated by 
Ui,... , w, Vi+i,. ..,Vjc. The vector w may lie outside of C, just as in the example in Figure 16 where 
w - (0,1). If If lies outside of C, then some of the Ci will have a negative sign in the decomposition 
(depending on whether the determinants detC and detC; have the same sign or not). Also, some of the 
Ci will be half-open. Which faces of the Ci are open can be determined via an entirely combinatorial 
case-analysis, as given in [42]. The algorithm then proceeds by recursion, constructing a decomposition 
of each of the Q in turn. 

The parameter that guarantees termination of this recurrence is the number of lattice points in the 
fundamental parallelepiped of C, which we call the index ind(C) of C. In a each recursive step, a given 
cone C is split into d new cones Q, until unimodular cones are obtained. This gives rise to a d-ary 
recursion tree T. The key property is that w is “short” enough so that the index of all d new cones is 
significantly smaller than the index of C. In fact, the index decreases so rapidly that the depth £ of T 
is guaranteed to be doubly logarithmic in the index of C, more precisely £ e 0{d ■ ln(ln(ind(C)))). The 
double logarithm is crucial, since ln(ind(C)) is polynomial in the encoding size of C. The number of 
terms in the Anal decomposition, i.e., the number of leaves of T, is then bounded by ln(ind(C))®^*'^*"^'^'^ 
If d is fixed, the number of terms is thus bounded by a polynomial in the encoding size of C. See [13, 
p. 140-141] for details. 

In summary, Barvinok’s algorithm provides us with the following result. 


Theorem 3.5 (Barvinok [16], see also [14, 33, 42]). Let C be a simplicial symbolic cone in dimension d 
with encoding size s. Then there exists an N e Z^o. unimodular d-dimensional half-open symbolic cones 
Cl,..., Cjv with the same apex as C and signs ai,...,ajvE{“lil} such that 

N 

[C] = J^ailCi] 

and, for any fixed d, N is bounded by a polynomial in s. Moreover, let r’jq,..., Vi^a denote the generators 
of the Ci and let Ui denote the unique lattice point in the fundamental parallelepiped ofCi. Then 


N 


= Y.ai 

i=l 


(1 - ■ ..." (1 - z’^td) 


and, again, the encoding size of this rational function is bounded by a polynomial in s, provided that d is 
fixed. 


Comparing Theorem 3.4 with Theorem 3.5 we And that Barvinok’s decomposition has the distinct 
advantage that the size of the rational function expressions it produces can be exponentially smaller 
than the expression obtained through explicit fundamental parallelepiped enumeration. This differ¬ 
ence can be particularly pronounced when the number of variables of the linear Diophantine system is 
small, but the coefflcients in the system are very large and give rise to cones in the decomposition that 
have large indices. (However, in Section 5 we are going to see an example of a class of linear Diophan¬ 
tine systems with large coefflcients where the indices of the cones produced in the previous step are 
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Still small.) The drawbacks of Barvinok decompositions are that they are determined by a non-trivial 
recursive algorithm, as opposed to the explicit formula of Theorem 3.4, and that the number of distinct 
factors (1 - z'') appearing in the denominators of the final rational function expression is much larger 
using the Barvinok approach than using fundamental parallelepiped enumeration. 

It is important to remark on the dimension d that drives the exponential growth of Barvinok’s de¬ 
composition. In the very first step of our algorithm, we apply the MacMahon lifting and transform a 
linear system with d variables and m constraints into a d-dimensional cone in {d + m)-dimensional 
space. However, the m extra dimensions are all eliminated during the second phase of the algorithm! 
Therefore the cones C for which we apply Barvinok’s algorithm all have dimension (at most) d and lie 
in d-dimensional space again. Thus, the MacMahon-style iterative elimination approach we employ 
does not lead to an exponential growth of the Barvinok decomposition, because we work with symbolic 
cones throughout the elimination and convert to rational function only after all extra variables have 
been eliminated. 

In practice, both approaches should be used in conjunction, as is also done in software packages 
like LattE [33]. If a symbolic cone C has a small index, fundamental parallelepiped enumeration is 
applied to directly obtain the rational function. If the index of C is large, Barvinok’s algorithm is used 
to recursively decompose C. However, instead of stopping the recurrence when ind(Q) = 1, i.e., when 
the cones C* become unimodular, the recurrence is stopped when the indices become small enough 
so that fundamental parallelepipeds can be explicitly enumerated efficiently. The threshold when the 
index becomes small enough is implementation dependent. 

3.4. Rational Function in Normal Form. No matter which mechanism for conversion to rational func¬ 
tions is used, the result is a rational function expression for <l>p, not a rational function in normal form. 
In a majority of applications, the desired output will be a rational function expression, as opposed to a 
rational function in normal form. In particular, if the linear Diophantine system has finitely many solu¬ 
tions, the normal form of the rational function will be a polynomial that explicitly lists all solutions, in 
which case it is often preferable to stick to a more compact expression. Even if the solution set is infinite, 
bringing the rational function in normal form - by bringing all summands on a common denominator, 
expanding and canceling terms - will typically destroy a lot of the structure of the rational function 
expression and may vastly increase its size. In fact, it is easy to construct families of cases where the 
size of the normal form is exponential in the size of the output produced by either of the conversion 
mechanisms presented above. 

That said, there are certainly cases where computing the normal form can involve helpful cancel¬ 
lations that indeed simplify the rational function expressions. In cases where this behavior can be 
expected, it may be useful to bring the rational function in normal form or to at least perform some 
simplifications. Normal form algorithms are widely available in current computer algebra systems and 
can readily be used for this purpose. We will therefore discuss this subject no further. We merely remark 
that in our implementation using the Sage computer algebra system, using Sage’s built-in mechanism 
for bringing the rational function in normal form often takes orders of magnitude longer than all the 
previous steps of the algorithm combined, due to the size of the expressions involved. 

4. Computational Complexity 

In this section, we give a detailed complexity analysis of Polyhedral Omega, providing bounds on the 
running time and space requirements of the algorithm. Polyhedral Omega is the first partition analysis 
algorithm for which this has been done. A key reason for this is that the geometric point of view enables 
us to prove stronger invariants than a naive complexity analysis would allow. In particular, we prove a 
structure theorem that implies strong bounds on the number of symbolic cones that can appear during 
elimination and on their encoding size. 

4.1. Structure Theorem. In order to prove a structure theorem for the symbolic cones that can appear 
during elimination inductively, we need to set up a suitable invariant that is preserved by el i m LastCoord. 

For brevity, we introduce some notation for matrix manipulations. Let Idix/ denote the rectangular 
(; X j)-matrix with ones on the diagonal and zeros everywhere else. Let /) £ 

denote the submatrix of a matrix A = {a^,v) n=i,...,m]v=i,...,n consisting of rows and columns 
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i.e., for and v = 1 ,..., /-A;+1 the entry 5^,v of the suhmatrix satisfies 

If b e is a vector, we abbreviate b^i^j) - i). In what follows we will assume a dense 

representation of matrices and vectors in order to simplify the presentation. We note, though, that in 
practice most of the matrices are sparse. We also assume constant time for accessing any element of a 
matrix. Generally, the estimates we obtain in this section are coarse. In particular, we work with naive 
bounds on the complexity of integer and matrix multiplication. A finer analysis would thus yield tighter 
bounds. 

Any simplicial cone can be described both in terms of its generators, as we have done in our sym¬ 
bolic cone representation, and in terms of a system of equations and inequalities. For the purpose of 
characterizing the cones that can appear during one run of elim, it will be convenient to work with 
inequality descriptions of symbolic cones. Given matrices A' e ■z^xid+m)^ ^ jdx{d+m) vectors 
b' E Z"* b" E such that 

= |x E 05^*+"* I A'x = h' and h"| 

is a d-dimensional simplicial cone in there is a unique set of d primitive integer vectors that 

generate K. Let V denote the matrix with these generators as columns, let e Q” denote the unique 
solution of A' q - b' and A''q = b” and let C = (V, q,0) denote the corresponding closed symbolic cone. 
Then K - [C] and we define C as the cone determined by the system A'x - b', b". Going one step 

further, we define flip(C) as the forward cone determined by the system A'x - b',A"x S b". For scalars 
a, f and vectors u, v we define 

and u>° v ^ Ui Vi for all i. 

Note that a superscript 1 reverses both the direction and the strictness of the inequality. Given this 
notation, we can observe that for flip(C) = (V, q, d) we have 

[flip(C)] = |x E I A'x = b' and A”x>° h"| 

where, without loss of generality, we assume that the columns vt of V are indexed such that the gener¬ 
ator Vi lies on the extreme ray of the cone that is given by letting all the constraints A”x'»° b" hold at 
equality, except for the i-th one. 

So, the openness vector o describes both, which facets of the cone are open and, correspondingly, 
which generators have been reversed by flip so that the cone points forward. One important observation 
is that no matter how a cone has been flipped in the past, the forward cone is uniquely determined by 
the original system. More precisely, 

(14) flip({x I A'x= b',A''x^° b”}) =flip({x | A'x - b',A”x^ b”}) 

for every vector o. 

With these preparations we can characterize all cones that can appear throughout the elimination 
process. 

Theorem 4.1. Let C - (V, q, o) be a symbolic cone with V e q e and o = 0 e {0,1}'* such 

that the columns of V are primitive and linearly independent and such that C is forward. Moreover, letC 
have an inequality description of the form 

[C] = |xElR'^+"' I (A -ldm)x=band{Id 0)xso| 

for some A E and b e Z™. Then, any symbolic cone C' = {V, q', o') produced at iteration i-Q,...,m 

o/elim(C) is a forward cone determined by a system of the form 

[C'l =flip(|xElR''+'""' I A'x^b',A"x^b"j) = |xe[R''+"'"' | A!x^ h',A!'x^°' b"j 

where A' = (A -Idm)(i_^_;)x(i_rf+m-/) ^ b' = h(i-m-r) e andA"x>b" isasubsys- 

tem consisting ofd rows of 

ftldx{d+m-i) ]x>[ 

(a .d+m-iV 

where e IR'* denotes the zero vector. In particular. A" e and b" e 
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The above theorem can be generalized to input cones C that have open faces, but this is not neces¬ 
sary for our purposes. 

Proof. The proof proceeds by induction over i. For i = 0, i.e., before elimLastCoord is applied for the first 
time, the claim is trivially true: The claim simply asserts that o) x S 0, (a -idm) x = h is an inequality 
description of C which is true by assumption. 

So we assume that the claim holds for some i and show that it then holds for / -i-1. Let C' - [V', q', o') 
be a d-dimensional cone produced at iteration i and let 

[C'l = jxElR''-""'"* I A!x^b'b"j 

be the system determining C', where Jf, b' ,A", b” are as described in the statement of the theorem and 
the columns of V are indexed to fit the order of the rows of A!' as described above. 

At iteration i -i-1 we intersect C' with the half-space given by xa+m-i ^ 0- By assumption on A',A", 
the variable xa+m-i has a non-zero coefficient only in the last row of Al and can therefore be seen as a 


slack variable. We rewrite 



^d+m-i ^ 0 

A 

Cl\X\ + ... + Cl^Xd Xd+ffi-i — 


A 

aiXi^...^adXd^b'j^_ 


where the a, are the first d entries of the last row of A'. Projecting away the coordinate x^+m-it we find 
that the polyhedron P - 7r([C'] n c jg given by the inequality description 

P = {x e I A'x = b', A"x b"j 

where A' = = (a b' - = h(i^m-!-i) and 

A"x b” is the system A"x b” with the additional constraint aiXi +... + a^ixa'^ b'^-i added. 

Consider first cases A and B from Section 3.2. As we observed, the cones produced by elimLastCoord 
are forward flips of the vertex cones of P. The vertex cones of a d-dimensional polyhedron in (d -i- m - 
i -1) -dimensional space are given by a system of m -1 -1 equations and d inequalities. In our case, the 
m - i- \ equations are precisely A'x = b' and the d inequalities are some subsystem of d inequalities 
chosen from A"x b". As shown in Section 3.2, case C is a deformation of case A. In particular, the 
cones produced in case C arise in the same way as in case A and hence are determined by the same type 
of system. 

We have now seen that any cone C produced at iteration / -i-1 is of the form 

(15) [C] = flip(|xe[R^+"'"'"^ I A'x= S',A"x3:"h"}-) = flip(|xE[R''+'”"'“i | A'x = h', A"x h"}-) 

where A!'x^° b” is a subsystem consisting of d rows of A"x 3;“ b" and the second equality follows from 
(14). The system (15) is precisely of the form asserted in the theorem, which completes the induction 
proof □ 

In particular. Theorem 4.1 applies to the symbolic cones C given by the MacMahon lifting. 

Lemma 4.2. Let A e and b e Z'”, then the cone C - {V, q, 0) produced by macmahon (A, b) has the 
inequality description 

[C] = {xeK''+'"|(A -Ini]x^b,[Id 0)xso}. 

Proof Recall y = | and q = ( The elements of C have the form q + Va for a e IR^q. We compute 

(idd 0 ) lq+ Va) - a^O and (a -im){q+ Va) - h which completes the proof □ 

To use Theorem 4.1 for our complexity analysis, we need to bound the size of a symbolic cone in 
terms of the size of its inequality description. 

Lemma 4.3. LetC be a d-dimensional forward symbolic cone in with primitive generators that 

is determined by the system 

[C] =flip(|xElR'*+'”"' I A'x= h',A"x3^ h"|). 
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Let the encoding size of all the entries in A!, b', A", b” be bounded above by s. Then the encoding size ofC 
is bounded above by 2{d+l){d+m-i]^\ogid+m - i)s+d with the size of each entry of V and q bounded 
by 

Li := 2{d + m- i)log{d + m- i)s. 

While Li is a function of d, m, i, s, we will in the following write just Li whenever the values of d, m, s 
are understood from context. 

Proof Let C - {V,q,o). The columns of V are of the form prim(j;) where the vectors v are solutions to 
Av - ej for A = j e _ Similarly, q is the solution to Aq - b where b- |y,|. All the entries 

of the V and q are therefore quotients of determinants det(Aj;) and det(A), where the matrices A^ are 
obtained by replacing a column of A by ej or b, as appropriate. 

By Hadamard’s inequality, | det(A) I ^ O; 11 1 12 ^ Oj 11 <2; 11 1 where the aj are the columns of A. There¬ 

fore, the encoding size ofdet(A) is at most [d+m-i) \og{d+m-i)s. The same holds for the determinants 
det(Ajfc). The encoding size of the quotient is therefore at most 2{d+m - i)\og[d+ m- i]s, which 

is thus also an upper bound on the encoding size of the entries of the vectors v and q. 

Note that the vectors v are rational, while prim(j;) are always integer. In particular, primCy) may 
actually be longer than v. However, the entries of prim(f) are factors of the integers det(Afc), whence 
their encoding size is bounded above by [d+m - i)\og{d+ m - i)s. Therefore, 2{d+m- i)\og{d+ m - i)s 
gives a uniform bound on the encoding size of the entries of V and q and the total the encoding size of 
C is at most 2(d +l](d+ m- i]^\og(d + m- i]s+ d. □ 

Putting these results together we obtain bounds on the number and encoding size of the symbolic 
cones produced during elimination. These bounds are much tighter than a naive analysis of the al¬ 
gorithm would imply, due to the geometric insight used. In particular, the encoding sizes of all cones 
produced during the iterative procedure can be bounded in terms of the encoding size of the original 
system. Here we use that rational simplicial cones have a unique representation as a symbolic cone 
with primitive generators (which are stored in sorted order). 

Corollary 4.4. Let C be as in Theorem 4.1. Then the following hold: 

(1) The number of symbolic cones produced after iteration i -0,...,m o/elim(C) is at most 

(2) For each cone C' = {V',q',o') produced after iteration i o/elim(C), the encoding size of the entries 
ofV' and q' is bounded by Li - 2[d + m - i) log(d + m - i) s, where s is a bound on the encoding 
size of each entry in the inequality description ofC. 

(3) In particular, for C - macmahon(A, b), the above holds ifs is a bound on the encoding size of each 
entry in A, b. 

Proof (1) The cones produced at iteration i each correspond to an inequality description as given in 
Theorem 4.1. Since there are only ways to select d rows from a matrix with d + i rows, there can be 
at most distinct cones. 

(2) At step i, we know by Theorem 4.1 there exists an inequality description of C' such that the co¬ 
efficients in this inequality description are bounded by 5 . By Lemma 4.3 the unique primitive sym¬ 
bolic cone representation C' = [V, q', o') has the property that the entries of V', q' are bounded by 
2{d+ m- i)\og{d+ m- i)s. 

(3) Follows from the (2) and Lemma 4.2. □ 

4.2. Complexity of Lifting and Eiimination. We can now give a complexity analysis of the algorithms 
in Section 3.1 and 3.2. 

Lemma 4.5 (Complexity of MacMahon Lifting, Algorithm 1). On input A e b e Z"* with entries of 

encoding size at most s, macmahon runs in O [dms + d^) time using^[dms + d^) space. 

Proof The MacMahon lifting is trivial and in essence only requires the time and space to read the input 
and write down the result. □ 

Lemma 4.6 (Complexity of Flip, Algorithm 2). Given a symbolic cone C - {V,q,o) withV eZ”^^, q^Q" 
and o E {0,1}*^, computingf\\p{C] andsgn(C) takes timeO[kn) time and no extra space. 
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Proof. The flip of a symbolic cone requires for each generator to run through its entries until we find 
the first non-zero entry. If it is negative, the multiply the rest of the entries and the sign of the cone by 
-1 and perform a xor operation on o. This requires O {kn) time and no extra space. □ 

Lemma 4.7 (Complexity of Prim, Algorithm 5). Given a matrix V e with entries of encoding size at 

mosts, prim runs in timeO[kns^) andOikns) space. 

Proof. To make each column primitive, we compute the greatest common divisor of all the entries in 
the column and divide. For every column, the greatest common divisor computation and the division 
take time O [ns^) each. □ 

Lemma 4.8 (Complexity of Eliminate Last Coordinate, Algorithm 3). Let i - \,...,m. Let C - {V,q,o) 
denote a symbolic cone satisfying the conditions given in Theorem 4.1 for one of the cones C' produced at 
iteration i- \ o/elim when applied to an input system A, b in which each entry has encoding size bounded 
abovebys. In particular, C isad-dimensionalconein{d+m-i + Vj-dimensionalspace. Then, computing 
el i m LastCoord (C) requires at most time 

0 (d^(d + m - i)^log(d) log(d + m - i)^s^] 
and space O {d[d + m - i + l)(.d + l)Lj_i). 

Proof. By Theorem 4.1 and Corollary 4.4, we know that the encoding sizes of the entries in V and q are 
bounded above by . 

Computing el i m LastCoord (C), we first construct the index set / in O (d) time, represented as a binary 
vector using d bits. Then we construct the vectors wj for the appropriate j’s. For each entry of wj we 
perform 6 multiplications, 1 subtraction, 1 greatest common divisor computation and 2 exact divisions, 
where all the initial quantities had encoding length bounded by L;-i. The time complexity for the com¬ 
putation of each entry is O {T;_i) and the space needed is 6Li_i -i-1. Thus, for the computation of wj we 
need 0[{d+ m - i+ l)L^_^) time and {d+ m - i + 1) (6L,_i -i-1) space. 

Computing an entry of the matrix VTl amounts to taking a linear combination of two elements from 
V with coefficients from Tl. The time complexity for each entry is O (i;_ J and the encoding size of an 
entry of of the linear combination is bounded by 2L;_i -i-1. This is repeated for al\d(.d+m-i + l) entries 
of the matrix VTf resulting in time complexity 0(d(d + m- i + l)L?_j). 

For each of the columns of the new matrix, we forget the last coordinate, which costs nothing and 
does not affect the bounds in the encoding size. Then we use prim to make the columns of this trun¬ 
cated {d+m-i]x d-matrix primitive, which takes time at most 0(d(d + m - by Lemma 4.7. The 

resulting matrix is the generator matrix of a symbolic cone after the i-th elimination step, and thus by 
Theorem 4.1 and Corollary 4.4 the encoding size of its entries are bounded by Lj. Flipping the resulting 
cone forward takes an additional time of at most Q>{d{d+ m - i)) by Lemma 4.6. In total, computing 
one symbolic cone takes time 

O ((2d -I- 1) (d -I- m - I -I- + 2d+ m - i + 2) ^'0[d{d+ m - i + . 

In order to make comparison of cones faster in the next step, we sort the generators of the symbolic 
cone lexicographically. This requires dlog(d) comparisons of vectors of encoding length [d + m - i]Li. 

Finally, inserting the at most d resulting symbolic cones in a dictionary mapping a symbolic cone to 
its multiplicity costs 0(dlog(d)(d -i- l)(d-i- m - i)L;), since the comparison of two symbolic cones costs 
O ((d -I- l)(d-i- m - i)Li). 

The total complexity of eliminating the last coordinate is thus 

0(d^(d + m- i) -i- log(d)L;)) ^ 0(d^(d + m- /)log(d)L?_]^) 

^ 0(d^(d + m- /)^log(d)log(d-i- m- i)^s^) 


where, in the first step, we use L;_i S Li and, in the last step, we substitute Lt-i -2{d+ m - i) log(d -i- 
m-i)s. □ 
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Lemma 4.9 (Complexity of Eliminate Coordinates, Algorithm 4). Let C- (V,q,o) denote a symbolic cone 
satisfying the conditions given in Theorem 4.1. Then, computing el i m (C) requires at most time 


O 


^ ^d+m-l^ 


[d + m)^rf^m"^log(d)^log(d + mY 


and space <[Si^s[‘^^J") id + mYdlogid + m) . 


Proof. At elimination step / - 1, we obtain { ) symbolic cones. For each of them we compute a 

dictionary containing d symbolic cones using Algorithm 3 in + m - /)^log(rf) log(t/ + m - iYs^) 

time. Then we insert each cone appearing in these dictionaries into a new dictionary that will contain at 
most cones. Comparing two cones produced in elimination step i costs O {d{d + m - i)L{]. Using 
appropriate data structures, the complexity of inserting an element in a dictionary of length is 
log(('*^*)) eO(/log(rf)) comparisons. Thus, inserting a cone in the dictionary costs 0(/rflog(d)(d+m- 
i)Li). We have to insert in the dictionary the d['^^y^) symbolic cones obtained after eliminating the 
i-th coordinate in the symbolic cones obtained in the previous step. In total, the complexity of 

obtaining the dictionary at elimination step i is 


id^{d+ m- i)‘^log(d)log(d+ m- i)s 


The cost of computing el im LastCoord( C) for each of the symbolic cones is {d+m- i)^log(d) 

log(d + m - iYs^) from Lemma 4.8. Thus, elimination step i has a total cost of 



'd+i-Y 


' 

'd+i- 1 

o 


id^log{d]{d + m- i)Li 

) 

^ 0 



O 


^ o 


^ d+i-l 

d+i-I 
d 


[id^{d+ m- irlog{d)logid+ m- i)s+ d^{d+ m- /)^log(d)log(t/+ m- irsy 


\ \ 


id {d + m- /) log(rf) log(d + m- i) s 


Iterating for i - 1,2,. ..,m results in time complexity 

yd+m-l\ 


O 


d 


id + mY d^ m"^ logidY logid + mYs' 


2 ..2 


Regarding the space used we observe that in each iteration we need to store at most cones and 

that each cone has encoding size at most Oidid+ m)Lo], which yields the desired bound. □ 


4.3. Complexity of Rational Function Conversion. Next, we analyze how long it takes to explicitly enu¬ 
merate the lattice points in the fundamental parallelepiped of a symbolic cone C = iV, q, o), using the 
Smith normal form as described in enumFundPar (Algorithm 7). The total number of lattice points we 
have to enumerate is D = 1 det(U)| which can be exponential in the encoding length of V. As the size of 
the output can be exponential in the size of the input, it makes sense to give an output-sensitive analysis 
of the algorithm. In particular, the following lemma states the preprocessing time of enumFundPar, i.e., 
the time it takes to perform a one-off calculation at the beginning of the algorithm, and then the delay, 
which is the time it takes to generate the next lattice point in the list of lattice points in the fundamental 
parallelepiped. It turns out that the preprocessing time and the delay are polynomial in the input size 
and independent of the output size. 


Lemma 4.10 (Complexity of Enumerate Fundamental Parallelepiped). Let C - iV,q,o) denote a d- 
dimensional symbolic cone in such that the encoding size of each entry in V and q is bounded by 
s. Listing the D = |det(V)| lattice points in the fundamental parallelepiped ofC using enumFundPar 
as defined in Algorithm 7 takes preprocessing time O [d^ (log(r/) -i- s)^) and has delay O [d‘^ (log(rf) -i- s)^) 
resulting in a total running time of 

0(rf®(log(d) -I- s)^ -I- Dd'^ilogid) + s)^) . 

The space required during preprocessing and generation of individual lattice points is polynomial in s 
and d and independent of the total number D of generated lattice points. 
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Proof. The preprocessing phase of enumFundPar consists of computing the Smith normal form and 
then computing q, W~^S'q, and - Vq^’^^'^ + Sj^q. Since we assume that C is full-dimensional, 

p does not have to he computed, hut can he assumed to be p = 0. These operations involving q are 
relatively expensive as we need to work with rational vectors. Bounds on the running time required 
for computing the Smith normal form S, U~^ and W~^, as well as bounds on the encoding size of the 
entries in these matrices are taken from [58, 59]. In total, this gives a bound of O (d® (log(d) -i- s)^) on the 
running time of the preprocessing step. 

To generate the list of all lattice points, one lattice point at a time, we do not first generate xe P and 
then compute W~^S'x, as this would be inefficient. Instead, we observe that the lattice points in P can 
be listed in such a way that two successive entries in the list differ by a unit vector. Thus we can list the 
lattice points W~^S'P by starting at 0 and then adding or subtracting one column of W~^ S' in each step. 
After the modulus has been taken, multiplication with V is required, though. It is also important to note 
that the vectors V{W~^S'x mod" and - V+ sicq are both integer and that the final division by sjc 
is guaranteed to produce an integer vector. Since the final result is guaranteed to lie in the fundamental 
parallelepiped of C, the encoding size of every entry in every single output vector is bounded hy s+ d. 
In total, this gives a bound of O (t/^(log(d) -i- s)^) on the time taken to generate the next point in the list. 

While generating lattice points, only one point in W~^S'P has to be stored at a time. Therefore the 
space required is independent of the total number of points generated. □ 


To analyze the complexity of computing a Barvinok decomposition, we cite the following result due 
to Barvinok [16,17,14], see [13, 31] for details. 

Lemma 4.11 (Complexity of Barvinok Decomposition [17]). Let C - (V, q, o) he a d-dimensional sym¬ 
bolic cone in and let s be a bound on the encoding size of all entries in V andq. LetD-AeX.{V) denote 
the number of lattice points in the fundamental parallelepiped ofC. Then a list of the symbolic cones 
Cl,..., Cjv as given in Theorem 3.5 can be computed in time t{d, s) which is a polynomial in s for fixed d. 
Moreover, the number of cones satisfies 

logrf 

iV ^ d ■ (logD) = d ■ (log D) 5 ’ 

which is a polynomial in s for fixed d, sincelogD ^ ds-t |log(d). 

Unfortunately, an explicit formula for the bound t{d, s) is not available in the literature. The deriva¬ 
tion of such a bound, while straightforward in principle, is out of the scope of this article. 


Proof. The time bound, originally obtained in [17], is given in ]31, Theorem 7.1.10]. We use the bound 
on the depth of the decomposition tree from [31, Lemma 7.1.9] to get 


i°gi°gp 

N^dl ‘“ 83 ^ 


1 +- 


logd 


.^(^loglogDjlogg^ ^ d(logD]'°^T-l = d(logD)^'°g^ 

as an upper bound on the number of unimodular cones in the final decomposition. 


□ 


The complexity of both fundamental parallelepiped enumeration and Barvinok decomposition de¬ 
pend on the number D - \ det(U) | of lattice points in the fundamental parallelepiped of a given symbolic 
cone C = [V, q, o). 

Lemma 4.12. Let a e Z^q be an upper bound on the absolute value of the entries of the original linear 
Diophantine system Ax ^ b. Then for each cone C - {V,q,o) output by e\\m{macmahon{A,b)] thenumber 
of lattice points in the fundamental parallelepiped ofC is bounded above by 

|det(l/)l ^ {da/\ 

Proof. Let C = (V, q, o) be a symbolic cone output by elim(macmahon(A, b)). By Theorem 4.1 and the 
proof of Lemma 4.3, V is of the form V - prim(A“^) where prim(M) denotes the matrix obtained by 
applying prim to the columns of M and A is a d x d integer matrix with entries that are bounded in 
absolute value by a. The inverse A~^ has rational entries with common denominator det(A). Thus 
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det(j4) A Ms an integer matrix and each column of det( A) A Ms a positive multiple of the corresponding 
column of prim(A“'^). So 

|det(prim(A"^))| ^ |det(det(A)A"^)| = ldet(A)|'^ ■ |det(A“^)| = 1 detlAll^^'^ 

Again, using the Hadamard hound we obtain | det(A)| ^ {da}'^. Combining these inequalities we get 
I det(l^)| ^ [da)'^^ as desired. □ 


4.4. Summeuy of Complexity Analysis. The results of the above complexity analysis are summarized 
in the following theorem. 

Theorem 4.13. Let A e and b e Z"* such that the encoding size of all entries in A and b is bounded 
above by s. Then the different variants of the algorithm described in Section 3 have the following running 
times. 

(1) Computing a list of symbolic cones Ci,..., C„ with multiplicities cri,..., such that 

[{x \ Ax^ b,x>0}]^Y. “i[Ml 


U5i>igmacmahon andeWm takes time at most 
O 


21 

d-\- m-\ 

S \ 



[d + m)^rf^m^log(d)^log(d + m)^ 


This is a polynomial in the parameters, provided at least one ofd or m are fixed. The number of 
points ni in the fundamental parallelepiped of cone Ci is bounded above by ni ^ {d2^)‘^ and the 
total number of cones is bounded above byn^ rf”*) ■ 

(2) After the list of symbolic cones has been computed, there are two ways to obtain a rational func¬ 
tion representation of the desired generating function | Ax»b,x»o} (-z) ■ 

(a) Using fundamental parallelepiped enumeration, a rational function expression of the form 


$ 


{x I Ax>b,x^0} 








i=i 


can be emitted iteratively as follows: The at and vij are given directly by the list of symbolic 
cones. After preprocessing time 0[nd^^\og{d)^s^] the at mostY. j n j vectorsuij canbelisted 
one at a time with delay O [d^ log{d)^ s^). In total, the time taken for computing the entire 
expression is thus at most 


O 


d-i-m' 

V ^ y 




(b) Using the Barvinok decomposition, a rational function expression of the form 


N 


0) 


{x I Ax^b,x>0} 






can be computed in time at most 


O 


d-\-m 

d 


t{d,2d\og{d)s) 


where t{d, s] is the expression from Lemma 4.11 which is a polynomial in s ifd is fixed. The 
number N of summands is at most 

d 


N ^ 


V ^ I 


<i(d2log(rf)s)<^-l'2)logd_ 


All the bounds stated in this theorem are rather coarse and could be improved by a refined analysis. 
Note also that, by construction, N ^ n, i.e., the Barvinok decomposition will produce always at least as 
many summands in the outer sum. Yet, typically N will be exponentially shorter than X; tt/. In (2a), the 
exponents Vj j will be the generators of the symbolic cones produced in (1) and thus edge directions 
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Elim 

RF 

Gn 

P 

G2 

Elliott/MacMahon 

R 



R 

R 

Omega2 

R 



R 

L 

Ell 

R 



R 

L 

CTEuclid 

R 


R 



PolyOmega-FP 

R 

L 




PolyOmega-BD 

R 

R 





Table 1. Iteration in partition analysis algorithms. See text for abbreviations. 


of the arrangement of hyperplanes given by the original system. In (2b), the set of exponents Wij will 
typically be much larger, since the wij arise from the Vij- through a recursive arithmetic procedure. 

Proof. (1) follows from Corollary 4.4 and Lemmas 4.5, 4.9, 4.12. For (2a), we combine Corollary 4.4 and 
Lemma 4.10. For (2b), we combine Lemmas 4.11 and 4.12. □ 


We conclude this complexity analysis with a couple of final remarks. (1) runs in polynomial time if 
either d or m is fixed. (2b) runs in polynomial time if d is fixed. (2a) does not run in polynomial time 
if d is fixed, due to the term 2* which arises from enumerating fundamental parallelepipeds which may 
be exponentially large. This means that using (1) and (2b) we obtain an algorithm that solves the rfsLDS 
problem in polynomial time, if the dimension of the problem is fixed. Thus, the Polyhedral Omega 
algorithm lies in the same complexity class as the best known polyhedral algorithms for solving rfsLDS. 
None of the previous partition analysis algorithms is known to have this property, see Section 5.2. 

A key bottleneck is the number of symbolic cones output by elim for which we have the upper bound 
upper bound is polynomial if either dor m are fixed, but can be exponential otherwise, e.g., 
=: A bottleneck of this type is inherent in all algorithms that have to touch each vertex of the 

polyhedron defined by the inequality system at least once (e.g., when using Brion decompositions): By 
McMullen’s Upper Bound Theorem [51], we know that the maximal number of vertices that a polytope 
with a given number of facets can have is realized by duals of cyclic polytopes. In our context, this yields 
inequality systems with 


'd+ m- 

I - I 

' [2 J 


r 1 \ 

2 


(d + m-l- I'^ll 

' ' 

1 




1 d-1 1 

) 


^ J 


vertices [66]. It is important to remark that the actual number of symbolic cones is typically much lower 
than ( but can be larger than the number of vertices, even in the case of simple polyhedra. See 

Section 5.3 for further discussion and open research questions related to this phenomenon. 


5. Comparison with Related Work 

In this section, we compare Polyhedral Omega with previous algorithms from both partition analysis 
and polyhedral geometry. 

5.1. Speedup via Symbolic Cones. Partition analysis algorithms are characterized by the iterative ap¬ 
plication of a collection of explicitly given formulas. However, iteration can appear in the various algo¬ 
rithms on many different levels and in many different forms, which can have important performance 
implications. 

Table 1 gives an overview of the iterative structure of the classic algorithm by Elliott/MacMahon, 
the Omega2 algorithm by Andrews-Paule-Riese [5], Xin’s algorithms Ell [63] and CTEuclid [65] and 
the two variants of our Polyhedral Omega algorithm, using fundamental parallelepiped enumeration 
(PolyOmega-FP) and Barvinok decompositions (PolyOmega-BD). Iteration can take either the form of 
a recurrence (R) or the form of a loop (L). By a loop we mean explicit formulas involving a long summa¬ 
tion. The defining feature of all partition analysis algorithms is the iterative elimination (Elim) of one 
variable at a time using a recursion that splits one rational function/cone into several rational func¬ 
tions/cones. In analogy with geometry, we will call the binomial factors in the denominator of a ra¬ 
tional function generators. Within the elimination of one variable, the algorithms have very different 
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Structure. Elliott’s algorithm uses an outer recurrence over pairs of generators (P) as well as an inner 
recurrence, similar to the Euclidean algorithm, to apply to one pair of generators (G2). A key con¬ 
tribution of OmegaZ was to replace the inner recurrence by an explicit formula (i.e., a loop) leading to 
a significant speedup. Xin’s Ell algorithm has a similar structure. In contrast, both Polyhedral Omega 
and Xin’s CTEuclid algorithm do not iterate over pairs of generators, but apply to all n generators 
at once (Gn). While CTEuclid uses a recursive procedure in this phase. Polyhedral Omega does not 
need to use any iteration at all here, due to the use of symbolic cones. In Polyhedral Omega, the cor¬ 
responding iteration can happen outside the elimination recurrence, when converting symbolic cones 
to rational functions (RE). The fact that in Polyhedral Omega it is not necessary to iterate at all inside 
a single elimination can already provide an exponential speedup in comparison with other partition 
analysis methods, as we will see in a concrete example below. The conversion to rational functions can 
be performed after elimination is complete using either fundamental parallelepiped enumeration (i.e., 
an explicit summation formula) or the Barvinok decomposition (i.e., a recursive procedure). 

Polyhedral Omega is the only partition analysis algorithm for which a comprehensive complexity 
analysis is available. Furthermore, Polyhedral Omega is the only partition analysis algorithm whose 
elimination phase runs in polynomial time if either the number of variables or the number of con¬ 
straints is fixed, and, if Barvinok decompositions are used, it is the only partition analysis algorithm 
that, in total, runs in polynomial time if the dimension is fixed. These are not just better bound due to a 
refined analysis that exploits the geometric structure, but genuine exponential speedups, even if Barvi¬ 
nok decompositions are not used. To illustrate these performance gains, we will discuss two structural 
differences between the algorithms that lead to exponential improvements on large classes of instances. 

First, it is generally desirable to handle all generators at once, instead of iterating over pairs of gen¬ 
erators. OmegaZ, for example, acts on a rational function by taking a product of two generators with 
positive last coordinate and splitting it into a sum of two rational functions that each have one genera¬ 
tor with positive last coordinate and one generator with zero last coordinate. Now, assume d is even and 
p is a rational function with d generators, half of which have positive and half of which have negative 
last coordinate. If we recursively apply the OmegaZ rule until each summand has at most one generator 
with positive last coordinate (which is a base case of the OmegaZ algorithm), we end up with a total of 
2 ^/ 2 -! summands. In contrast, Polyhedral Omega’s elimLastCoord produces at most d+ \ summands in 
the worst case. 

Second, and most importantly, using symbolic cones during elimination can lead to exponentially 
more compact intermediate representations, which in turn implies exponential speedups during the 
elimination phase, compared with algorithms working with rational functions. This can save a lot of 
work in total, as soon as several variables need to be eliminated, as the following example shows. Let <21 
be any large positive integer and let 02 he any positive integer coprime to . Let bi , £>2 denote positive 
integers such that <21 ^2 “ <22 ^1 = 1 • Now consider the linear Diophantine system S consisting of the two 
constraints bzXi -biXz^O and -a 2 X 2 + aiXi'^ 1. If we solve S iteratively by, in the first step, computing 
a full rational function representation of the set of all Diophantine solutions x ^ 0 to the inequality 
- 02 X 2 + aiXiS: 0, we will typically obtain a rational function expression such as 




or a lifted version thereof, which has ai terms in the numerator. This is because the fundamental par¬ 
allelepiped of the corresponding cone 1 )^ contains <21 lattice points. Using rational function 

representations, the first elimination has thus increased the encoding size of the problem exponen¬ 
tially, which will slow down subsequent eliminations correspondingly. Using symbolic cones, on the 
other hand, the encoding size does not increase during elimination: Polyhedral Omega starts with the 
symbolic cone 


Co = '^r(( I _t):0). 

\ —CI2 Cl\ ) 
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The first elimination produces 

I ai 0 \ f ^ ' 

Ci='^r( 1 ; “1 ). 

U -fcJ -h. 

V aj / 

As expected, the projection of Q onto the first two coordinates is the solution set of the system x 0 
and -a 2 X 2 + aiXi 1. Note that, from Co to Ci, the determinant of the cone has increased from 1 to tti 
hut the encoding size has not. The second elimination then produces 

which is precisely the solution set of the system S. The cone C 2 has again determinant 1, so that the fun¬ 
damental parallelepiped can be enumerated in just a single iteration and the resulting rational function 
expression is short. Since we used symbolic cones we did not have to enumerate the fundamental par¬ 
allelepiped explicitly at the intermediate stage and can get to the short end-result quickly. Rational 
function representations, however, produce an exponential blowup at the intermediate stage that can 
prevent partition analysis algorithms from reaching the short final result. 

Of course there are shorter expressions for the rational function p given in (16), which can be ob¬ 
tained, e.g., using Barvinok decomposition or methods from [22, 65]. This would avoid an exponential 
blowup in one elimination step. But the increase in representation size in a single elimination step 
would still be larger than using symbolic cones. It is thus preferable to apply Barvinok techniques only 
after elimination has been completed, as we do in Polyhedral Omega. Initial experiments with our 
implementation of Polyhedral Omega [23] confirm these benefits of symbolic cone representations in 
practice. 

As the above example shows, symbolic cones are the key for obtaining Polyhedral Omega’s improve¬ 
ments in running time. However, symbolic cones are not tied to the particulars of our algorithm. On the 
contrary, symbolic cones have great potential for further applications in partition analysis and beyond. 
One particularly promising direction for future research is to develop a symbolic cone version of partial 
fraction decomposition, following the polyhedral interpretation of PFD given in Section 2.6. This might 
lead in particular to a polyhedral version of Xin’s CTEuclid algorithm with improved performance when 
several variables are eliminated. 

5.2. First Polynomial-Time Algorithm from Partition Analysis Family. Polyhedral Omega is the first 
algorithm from the partition analysis family for which a comprehensive complexity analysis is available 
(Section 4). Geometric insight played a crucial role in obtaining strong bound on the running time. 
Most importantly, by Theorem 4.13, Polyhedral Omega runs in polynomial time if the dimension is 
fixed, provided that the Barvinok decomposition is used. It is the first partition analysis algorithm with 
this property. 

This latter statement requires some discussion. First, it is important to draw the following careful 
distinction: The Barvinok decomposition is an algorithm for computing a short signed decomposition 
of a simplicial cone into unimodular simplicial cones, which runs in polynomial time if the dimension is 
fixed. More generally, however, there is Barvinok’s algorithm for solving the rfsLDS problem. Barvinok’s 
algorithm makes crucial use of Barvinok decompositions, but performs a whole range of other tasks as 
well, including but not limited to computing vertices and edge directions as well as triangulation. 

As we mentioned in Section 2.5, [65, Section 3] describes a procedure for reducing the problem of ap¬ 
plying to a rational function to the rfsLDS problem. The rfsLDS problem itself, though, is solved by 
using Barvinok’s algorithm as a black box. In our view, the procedure from [65, Section 3] does therefore 
not fall into the domain of partition analysis algorithms. In contrast to [65, Section 3], Polyhedral Omega 
does not use all of Barvinok’s algorithm as a black box, but only Barvinok decompositions. The Barvi¬ 
nok decomposition itself can be described in terms of a recursive formula [42], very much in the spirit 
of partition analysis. Thus, we argue. Polyhedral Omega with Barvinok decompositions (PolyOmega- 
BD) is the first rfsLDS algorithm that follows the partition analysis paradigm and achieves polynomial 
running time in fixed dimension. 

The CTEuclid algorithm given in [65, Section 4] is quite separate from the algorithm from [65, Section 
3] mentioned above, and does not use Barvinok decompositions at all. Xin raises the question whether 
CTEuclid runs in polynomial time if just a single variable is eliminated. This question remains open. 
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Figure 17. This example shows how zero sums of symbolic cones can arise dur¬ 
ing elimination. Consider the area S above and including the “W”-shaped line. (A) 
gives a symbolic cone decomposition [S] = [Ci] -i- [C 2 ] - [C 3 ]. Intersecting each of 
these symbolic cones with the half-space above the horizontal axis gives the de¬ 
composition [S n H^] = [Ci,i] - [Cl, 2 ] + [C 2 ,i] - [C 2 , 2 ] - [C 3 ] shown in (B). Note that 
[0] = [C 2 ,i] - [C 3 ] - [Cl, 2 ], but these cones do not cancel when collecting terms. Espe¬ 
cially in higher dimension, situations like these can arise even if the starting set S is 
convex and no cut goes through the apex of a cone. 


However, as Xin acknowledges, even if the answer was positive, polynomiality would likely break down 
as soon as several variables are eliminated one after another. Here symbolic cones could be of help as 
explained in Section 5.1. 

5.3. Comparison with Polyhedral Methods. In comparison with polyhedral methods, the main benefit 
of Polyhedral Omega is simplicity. Polyhedral algorithms for rfsLDS typically make use of a number of 
non-trivial and highly-specialized algorithms for enumerating vertices and edges of a polyhedron and 
for triangulating a cone, since these tasks are essential for decomposing a polyhedron into simplicial 
cones. In contrast. Polyhedral Omega obtains such a decomposition implicitly through the iterative 
application of a few simple rules. This fact is the main insight transferred from partition analysis to 
polyhedral geometry as part of this research project. Crucially, this simple iterative approach still results 
in an algorithm that lies in the same complexity class as the best knovm polyhedral methods for solving 
rfsLDS: Polyhedral Omega runs in polynomial time in fixed dimension. 

However, there is one interesting regard in which there is a conceptual gap between the performance 
of Polyhedral Omega and polyhedral algorithms using specialized methods for vertex and edge enumer¬ 
ation. Applying our iterative approach, we obtain a representation 

N 

[P] = ^a;[Q] 

i=i 

of a polyhedron P as a sum of (indicator functions of) symbolic cones Cjwith multiplicities ai e Z. As 
we discussed in Section 3.2 normalizing and collecting terms [C,] is crucial for performance. However, 
over the course of the algorithm, there may still arise groups of symbolic cones that sum to the empty 
set, i.e., 

l 

ai[Ci\ = [0], 

i=k 

without our algorithm being able to recognize this fact and simplify the sum accordingly. To see how 
such zero sums can arise, consider the example given in Figure 17. The underlying issue here is that 
on the rational function level such zero sums could be recognized by bringing the rational function 
expression in normal form. This calls for the development of a symbolic method for bringing symbolic 
cone expressions into a normal form. Such a method could speed up Polyhedral Omega significantly 
and is a promising topic for future research. 

6. Summary AND Future Research 

In this article, we have given a detailed polyhedral interpretation of existing partition analysis algo¬ 
rithms solving the rfsLDS problem (Section 2) and used these insights to develop Polyhedral Omega 
(Section 3), a new algorithm solving rfsLDS which combines the polyhedral and the partition analysis 
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approaches in a common framework. A key insight is the use of symbolic cones instead of rational func¬ 
tions, which allows us to retain the simple iterative approach based on explicit formulas that is inher¬ 
ent to partition analysis algorithms and still achieve a running time that had previously been attained 
only by polyhedral algorithms (Sections 3.2 and 5.1). For conversion from symbolic cones to rational 
functions, we present two methods (Section 3.3). The first, based on fundamental parallelepiped enu¬ 
meration (Theorem 3.4), has the advantage of being described by a simple closed formula, while the 
second, based on Barvinok decompositions (Theorem 3.5), is a recursive algorithm that yields short 
rational function expressions. When using Barvinok decomposition. Polyhedral Omega runs in poly¬ 
nomial time in fixed dimension (Theorem 4.13). Thus it lies in the same complexity class as the fastest 
known polyhedral algorithms for rfsLDS and it is the first algorithm from the partition analysis family to 
do so (Section 5.2). Moreover the geometric point of view allows us to do the first comprehensive com¬ 
plexity analysis of any partition analysis algorithm, and obtain strong bounds on the running time from 
the geometric invariants we establish (Theorem 4.1). Compared to previous polyhedral methods, the 
key advantage of Polyhedral Omega is its simplicity: instead of using several sophisticated algorithmic 
tools from polyhedral geometry, a simple recursive algorithm based on a few explicit formulas suffices. 

Our results point to several promising directions for future research. 

On the practical side, the next task is to complete our implementation of Polyhedral Omega [23]. 
So far, we have implemented the elimination algorithm using symbolic cones as well as rational func¬ 
tion conversion using fundamental parallelepiped enumeration on the Sage system. The next step is 
the implementation of Barvinok decompositions on this platform. This would enable a comprehen¬ 
sive benchmark comparison of the different rfsLDS algorithms available today. In this context, it would 
also be desirable to do a detailed complexity analysis of the Barvinok decomposition, giving an explicit 
formula for an upper bound on the running time of that algorithm. While such an analysis is straight¬ 
forward in principle, no such formula has been published so far. 

On the theoretical side, the most interesting direction for future research is a detailed study of sym¬ 
bolic cones. Symbolic cones form an interesting algebraic structure closely related to, but distinct from 
both, the rational function field and the algebra of polyhedra [13]. In particular, a good definition of a 
normal form of a symbolic cone expression along with an efficient algorithm for bringing a symbolic 
cone expression into normal form would be of great use, both in the context of Polyhedral Omega (see 
Section 5.1) as well as for working with rational functions in general. Moreover, symbolic cones may 
be of use outside the context of Brion/Lawrence-Varchenko decompositions. Especially our polyhedral 
interpretation of partial fraction decomposition in Section 2.6 calls for closer investigation. A symbolic 
cone variant of partial fraction decomposition could, for example, lead to a polyhedral version of the 
CTEuclid algorithm. 

Finally, a key feature of partition analysis algorithms has always been that they lend themselves very 
well for use in induction proofs, due to their being based on the recursive application of explicit rules. 
Polyhedral Omega promises to be of great use in both manual and automatic induction proofs, since the 
intermediate results obtained are given as symbolic cones and therefore more compact and easy to ma¬ 
nipulate. We are especially looking forward to combining Polyhedral Omega with automatic induction 
provers, as this may lead to algorithms for proving structural results about infinite families of polytopes 
of varying dimension and even to algorithms for manipulating infinite-dimensional polyhedra. A first 
milestone in this direction would be to finally realize one of the visions that motivated Andrews-Paule- 
Riese to revive partition analysis in the first place: a fully automatic proof of the Lecture Hall Partition 
Theorem [4, 21] in full generality, with the dimension treated as a symbolic variable. 
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